text:
# Analog Out API

#------------------------------------ Configuration
function analogOutCount(hdwf::Hdwf)
  pcChannel = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutCount,libdwf),Cint,(Hdwf,Ptr{Cint}),
    hdwf,pcChannel)
  succ == 0 && error("Error calling FDwfAnalogOutCount.")
  return pcChannel[]
end

function analogOutMasterSet(hdwf::Hdwf,idxChannel::Int32,idxMaster::Int32)
  succ = ccall((:FDwfAnalogOutMasterSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,idxMaster)
  succ == 0 && error("Error calling FDwfAnalogOutMasterSet.")
  return nothing
end

function analogOutMasterGet(hdwf::Hdwf,idxChannel::Int32)
  pidxMaster = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutMasterGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pidxMaster)
  succ == 0 && error("Error calling FDwfAnalogOutMasterGet.")
  return pidxMaster[]
end

function analogOutTriggerSourceInfo(hdwf::Hdwf,idxChannel::Int32)
  pfstrigsrc = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutTriggerSourceInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pfstrigsrc)
  succ == 0 && error("Error calling FDwfAnalogOutTriggerSourceInfo.")
  return pfstrigsrc[]
end

function analogOutTriggerSourceSet(hdwf::Hdwf,idxChannel::Int32,trigsrc::TRIGSRC)
  succ = ccall((:FDwfAnalogOutTriggerSourceSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,Cint(trigsrc))
  succ == 0 && error("Error calling FDwfAnalogOutTriggerSourceSet.")
  return nothing
end

function analogOutTriggerSourceGet(hdwf::Hdwf,idxChannel::Int32)
  ptrigsrc = Ref{Cuchar}()
  succ = ccall((:FDwfAnalogOutTriggerSourceGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cuchar}),
    hdwf,idxChannel,ptrigsrc)
  succ == 0 && error("Error calling FDwfAnalogOutTriggerSourceGet.")
  return TRIGSRC(ptrigsrc[])
end

function analogOutRunInfo(hdwf::Hdwf,idxChannel::Int32)
  psecMin = Ref{Cdouble}(); psecMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutRunInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,psecMin,psecMax)
  succ == 0 && error("Error calling FDwfAnalogOutRunInfo.")
  return (psecMin[],psecMax[])
end

function analogOutRunSet(hdwf::Hdwf,idxChannel::Int32,secRun::Float64)
  succ = ccall((:FDwfAnalogOutRunSet,libdwf),Cint,(Hdwf,Cint,Cdouble),
    hdwf,idxChannel,secRun)
  succ == 0 && error("Error calling FDwfAnalogOutRunSet.")
  return nothing
end

function analogOutRunGet(hdwf::Hdwf,idxChannel::Int32)
  psecRun = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutRunGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,psecRun)
  succ == 0 && error("Error calling FDwfAnalogOutRunGet.")
  return psecRun[]
end

function analogOutRunStatus(hdwf::Hdwf,idxChannel::Int32)
  psecRun = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutRunStatus,libdwf),Cint,(Hdwf,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,psecRun)
  succ == 0 && error("Error calling FDwfAnalogOutRunStatus.")
  return psecRun[]
end

function analogOutWaitInfo(hdwf::Hdwf,idxChannel::Int32)
  psecMin = Ref{Cdouble}(); psecMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutWaitInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,psecMin,psecMax)
  succ == 0 && error("Error calling FDwfAnalogOutWaitInfo.")
  return (psecMin[],psecMax[])
end

function analogOutWaitSet(hdwf::Hdwf,idxChannel::Int32,secWait::Float64)
  succ = ccall((:FDwfAnalogOutWaitSet,libdwf),Cint,(Hdwf,Cint,Cdouble),
    hdwf,idxChannel,secWait)
  succ == 0 && error("Error calling FDwfAnalogOutWaitSet.")
  return nothing
end

function analogOutWaitGet(hdwf::Hdwf,idxChannel::Int32)
  psecWait = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutWaitGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,psecWait)
  succ == 0 && error("Error calling FDwfAnalogOutWaitGet.")
  return psecWait[]
end

function analogOutRepeatInfo(hdwf::Hdwf,idxChannel::Int32)
  pnMin = Ref{Cint}(); pnMax = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutRepeatInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cint},Ptr{Cint}),
    hdwf,idxChannel,pnMin,pnMax)
  succ == 0 && error("Error calling FDwfAnalogOutRepeatInfo.")
  return (pnMin[],pnMax[])
end

function analogOutRepeatSet(hdwf::Hdwf,idxChannel::Int32,cRepeat::Int32)
  succ = ccall((:FDwfAnalogOutRepeatSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,cRepeat)
  succ == 0 && error("Error calling FDwfAnalogOutRepeatSet.")
  return nothing
end

function analogOutRepeatGet(hdwf::Hdwf,idxChannel::Int32)
  pcRepeat = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutRepeatGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pcRepeat)
  succ == 0 && error("Error calling FDwfAnalogOutRepeatGet.")
  return pcRepeat[]
end

function analogOutRepeatStatus(hdwf::Hdwf,idxChannel::Int32)
  pcRepeat = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutRepeatStatus,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pcRepeat)
  succ == 0 && error("Error calling FDwfAnalogOutRepeatStatus.")
  return pcRepeat[]
end

function analogOutRepeatTriggerSet(hdwf::Hdwf,idxChannel::Int32,fRepeatTrigger::Bool)
  succ = ccall((:FDwfAnalogOutRepeatTriggerSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,Cint(fRepeatTrigger))
  succ == 0 && error("Error calling FDwfAnalogOutRepeatTriggerSet.")
  return nothing
end

function analogOutRepeatTriggerGet(hdwf::Hdwf,idxChannel::Int32)
  pfRepeatTrigger = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutRepeatTriggerGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pfRepeatTrigger)
  succ == 0 && error("Error calling FDwfAnalogOutRepeatTriggerGet.")
  return Bool(pfRepeatTrigger[])
end

function analogOutModeSet(hdwf::Hdwf,idxChannel::Int32,mode::DWFANALOGOUTMODE)
  succ = ccall((:FDwfAnalogOutModeSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,Cint(mode))
  succ == 0 && error("Error calling FDwfAnalogOutModeSet.")
  return nothing
end

function analogOutModeGet(hdwf::Hdwf,idxChannel::Int32)
  pmode = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutModeGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pmode)
  succ == 0 && error("Error calling FDwfAnalogOutModeGet.")
  return DWFANALOGOUTMODE(pmode[])
end

function analogOutIdleInfo(hdwf::Hdwf,idxChannel::Int32)
  pfsidle = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutIdleInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pfsidle)
  succ == 0 && error("Error calling FDwfAnalogOutIdleInfo.")
  return pfsidle[]
end

function analogOutIdleSet(hdwf::Hdwf,idxChannel::Int32,idle::DWFANALOGOUTIDLE)
  succ = ccall((:FDwfAnalogOutIdleSet,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,Cint(idle))
  succ == 0 && error("Error calling FDwfAnalogOutIdleSet.")
  return nothing
end

function analogOutIdleGet(hdwf::Hdwf,idxChannel::Int32)
  pidle = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutIdleGet,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pidle)
  succ == 0 && error("Error calling FDwfAnalogOutIdleGet.")
  return DWFANALOGOUTIDLE(pidle[])
end

function analogOutNodeInfo(hdwf::Hdwf,idxChannel::Int32)
  pfsnode = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutNodeInfo,libdwf),Cint,(Hdwf,Cint,Ptr{Cint}),
    hdwf,idxChannel,pfsnode)
  succ == 0 && error("Error calling FDwfAnalogOutNodeInfo.")
  return pfsnode[]
end

function analogOutNodeEnableSet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,fEnable::Bool)
  succ = ccall((:FDwfAnalogOutNodeEnableSet,libdwf),Cint,(Hdwf,Cint,Cint,Cint),
    hdwf,idxChannel,Cint(node),Cint(fEnable))
  succ == 0 && error("Error calling FDwfAnalogOutNodeEnableSet.")
  return nothing
end

function analogOutNodeEnableGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pfEnable = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutNodeEnableGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cint}),
    hdwf,idxChannel,Cint(node),pfEnable)
  succ == 0 && error("Error calling FDwfAnalogOutNodeEnableGet.")
  return Bool(pfEnable[])
end

function analogOutNodeFunctionInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pfsfunc = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutNodeFunctionInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cint}),
    hdwf,idxChannel,Cint(node),pfsfunc)
  succ == 0 && error("Error calling FDwfAnalogOutNodeFunctionInfo.")
  return pfsfunc[]
end

function analogOutNodeFunctionSet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,func::FUNC)
  succ = ccall((:FDwfAnalogOutNodeFunctionSet,libdwf),Cint,(Hdwf,Cint,Cint,Cuchar),
    hdwf,idxChannel,Cint(node),Cuchar(func))
  succ == 0 && error("Error calling FDwfAnalogOutNodeFunctionSet.")
  return nothing
end

function analogOutNodeFunctionGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pfunc = Ref{Cuchar}()
  succ = ccall((:FDwfAnalogOutNodeFunctionGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cuchar}),
    hdwf,idxChannel,Cint(node),pfunc)
  succ == 0 && error("Error calling FDwfAnalogOutNodeFunctionGet.")
  return FUNC(pfunc[])
end

function analogOutNodeFrequencyInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  phzMin = Ref{Cdouble}(); phzMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeFrequencyInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,node,phzMin,phzMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodeFrequencyInfo.")
  return (phzMin[],phzMax[])
end

function analogOutNodeFrequencySet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,hzFrequency::Float64)
  succ = ccall((:FDwfAnalogOutNodeFrequencySet,libdwf),Cint,(Hdwf,Cint,Cint,Cdouble),
    hdwf,idxChannel,node,hzFrequency)
  succ == 0 && error("Error calling FDwfAnalogOutNodeFrequencySet.")
  return nothing
end

function analogOutNodeFrequencyGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  phzFrequency = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeFrequencyGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,node,phzFrequency)
  succ == 0 && error("Error calling FDwfAnalogOutNodeFrequencyGet.")
  return phzFrequency[]
end

function analogOutNodeAmplitudeInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pMin = Ref{Cdouble}(); pMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeAmplitudeInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,node,pMin,pMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodeAmplitudeInfo.")
  return (pMin[],pMax[])
end

function analogOutNodeAmplitudeSet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,vAmplitude::Float64)
  succ = ccall((:FDwfAnalogOutNodeAmplitudeSet,libdwf),Cint,(Hdwf,Cint,Cint,Cdouble),
    hdwf,idxChannel,node,vAmplitude)
  succ == 0 && error("Error calling FDwfAnalogOutNodeAmplitudeSet.")
  return nothing
end

function analogOutNodeAmplitudeGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pvAmplitude = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeAmplitudeGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,node,pvAmplitude)
  succ == 0 && error("Error calling FDwfAnalogOutNodeAmplitudeGet.")
  return pvAmplitude[]
end

function analogOutNodeOffsetInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pMin = Ref{Cdouble}(); pMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeOffsetInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,node,pMin,pMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodeOffsetInfo.")
  return (pMin[],pMax[])
end

function analogOutNodeOffsetSet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,vOffset::Float64)
  succ = ccall((:FDwfAnalogOutNodeOffsetSet,libdwf),Cint,(Hdwf,Cint,Cint,Cdouble),
    hdwf,idxChannel,node,vOffset)
  succ == 0 && error("Error calling FDwfAnalogOutNodeOffsetSet.")
  return nothing
end

function analogOutNodeOffsetGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pvOffset = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeOffsetGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,node,pvOffset)
  succ == 0 && error("Error calling FDwfAnalogOutNodeOffsetGet.")
  return pvOffset[]
end

function analogOutNodeSymmetryInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  ppercentageMin = Ref{Cdouble}(); ppercentageMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeSymmetryInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,node,ppercentageMin,ppercentageMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodeSymmetryInfo.")
  return (ppercentageMin[],ppercentageMax[])
end

function analogOutNodeSymmetrySet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,percentageSymmetry::Float64)
  succ = ccall((:FDwfAnalogOutNodeSymmetrySet,libdwf),Cint,(Hdwf,Cint,Cint,Cdouble),
    hdwf,idxChannel,node,percentageSymmetry)
  succ == 0 && error("Error calling FDwfAnalogOutNodeSymmetrySet.")
  return nothing
end

function analogOutNodeSymmetryGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  ppercentageSymmetry = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodeSymmetryGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,node,ppercentageSymmetry)
  succ == 0 && error("Error calling FDwfAnalogOutNodeSymmetryGet.")
  return ppercentageSymmetry[]
end

function analogOutNodePhaseInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pdegreeMin = Ref{Cdouble}(); pdegreeMax = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodePhaseInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),
    hdwf,idxChannel,node,pdegreeMin,pdegreeMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodePhaseInfo.")
  return (pdegreeMin[],pdegreeMax[])
end

function analogOutNodePhaseSet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,degreePhase::Float64)
  succ = ccall((:FDwfAnalogOutNodePhaseSet,libdwf),Cint,(Hdwf,Cint,Cint,Cdouble),
    hdwf,idxChannel,node,degreePhase)
  succ == 0 && error("Error calling FDwfAnalogOutNodePhaseSet.")
  return nothing
end

function analogOutNodePhaseGet(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pdegreePhase = Ref{Cdouble}()
  succ = ccall((:FDwfAnalogOutNodePhaseGet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble}),
    hdwf,idxChannel,node,pdegreePhase)
  succ == 0 && error("Error calling FDwfAnalogOutNodePhaseGet.")
  return pdegreePhase[]
end

function analogOutNodeDataInfo(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  pnSamplesMin = Ref{Cint}(); pnSamplesMax = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutNodeDataInfo,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cint},Ptr{Cint}),
    hdwf,idxChannel,node,pnSamplesMin,pnSamplesMax)
  succ == 0 && error("Error calling FDwfAnalogOutNodeDataInfo.")
  return (pnSamplesMin[],pnSamplesMax[])
end

function analogOutNodeDataSet(
    hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,rgdVoltData::Vector{Float64},len::Int32)
  succ = ccall((:FDwfAnalogOutNodeDataSet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Cint),
    hdwf,idxChannel,node, rgdVoltData,len)
  succ == 0 && error("Error calling FDwfAnalogOutNodeDataSet.")
  return rgdVoltData
end

function analogOutNodeDataSet(
    hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,rgdVoltData::Vector{Float64})
  succ = ccall((:FDwfAnalogOutNodeDataSet,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Cint),
    hdwf,idxChannel,node,rgdVoltData,Cint(length(rgdVoltData)))
  succ == 0 && error("Error calling FDwfAnalogOutNodeDataSet.")
  return rgdVoltData
end

#------------------------------------ Control
function analogOutReset(hdwf::Hdwf,idxChannel::Int32)
  succ = ccall((:FDwfAnalogOutReset,libdwf),Cint,(Hdwf,Cint),
    hdwf,idxChannel)
  succ == 0 && error("Error calling FDwfAnalogOutReset.")
  return nothing
end

function analogOutConfigure(hdwf::Hdwf,idxChannel::Int32,fStart::Bool)
  succ = ccall((:FDwfAnalogOutConfigure,libdwf),Cint,(Hdwf,Cint,Cint),
    hdwf,idxChannel,Cint(fStart))
  succ == 0 && error("Error calling FDwfAnalogOutConfigure.")
  return nothing
end

function analogOutStatus(hdwf::Hdwf,idxChannel::Int32)
  psts = Ref{Cuchar}()
  succ = ccall((:FDwfAnalogOutStatus,libdwf),Cint,(Hdwf,Cint,Ptr{Cuchar}),
    hdwf,idxChannel,psts)
  succ == 0 && error("Error calling FDwfAnalogOutStatus.")
  return DWFSTATE(psts[])
end

function analogOutNodePlayStatus(hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE)
  cdDataFree = Ref{Cint}(); cdDataLost = Ref{Cint}(); cdDataCorrupted = Ref{Cint}()
  succ = ccall((:FDwfAnalogOutNodePlayStatus,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint}),
    hdwf,idxChannel,node,cdDataFree,cdDataLost,cdDataCorrupted)
  succ == 0 && error("Error calling FDwfAnalogOutNodePlayStatus.")
  return (cdDataFree[],cdDataLost[],cdDataCorrupted[])
end

function analogOutNodePlayData(
    hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,rgdVoltData::Vector{Float64},len::Int32)
  succ = ccall((:FDwfAnalogOutNodePlayData,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Cint),
    hdwf,idxChannel,node,rgdVoltData,len)
  succ == 0 && error("Error calling FDwfAnalogOutNodePlayData.")
  return rgdVoltData
end

function analogOutNodePlayData(
    hdwf::Hdwf,idxChannel::Int32,node::ANALOGOUTNODE,rgdVoltData::Vector{Float64})
  succ = ccall((:FDwfAnalogOutNodePlayData,libdwf),Cint,(Hdwf,Cint,Cint,Ptr{Cdouble},Cint),
    hdwf,idxChannel,node, rgdVoltData,Cint(length(rgdVoltData)))
  succ == 0 && error("Error calling FDwfAnalogOutNodePlayData.")
  return rgdVoltData
end

#### TODO ####
# // EExplorer channel 3&4 current/voltage limitation
# DWFAPI BOOL FDwfAnalogOutLimitationInfo(HDWF hdwf, int idxChannel, double *pMin, double *pMax);
# DWFAPI BOOL FDwfAnalogOutLimitationSet(HDWF hdwf, int idxChannel, double limit);
# DWFAPI BOOL FDwfAnalogOutLimitationGet(HDWF hdwf, int idxChannel, double *plimit);
# // needed for EExplorer, don't care for ADiscovery
# DWFAPI BOOL FDwfAnalogOutCustomAMFMEnableSet(HDWF hdwf, int idxChannel, BOOL fEnable);
# DWFAPI BOOL FDwfAnalogOutCustomAMFMEnableGet(HDWF hdwf, int idxChannel, BOOL *pfEnable);


text:
using Caesar
using RoME, Distributions
import IncrementalInference: getSample

using TransformUtils

mutable struct ZJoint <: FunctorPairwise
  Zij::Distribution
end
function getSample(el::ZJoint, N=1)
  return (rand(el.Zij, N),)
end
function (el::ZJoint)(res, userdata, idx, meas, xi, xj)
  Xi = SE3(xi[1:3,idx],Euler(xi[4:6,idx]))
  Xj = SE3(xj[1:3,idx],Euler(xj[4:6,idx]))

    # h=meas[1][idx]
    # hip = SE3([0,0,0.5], Euler(0,0,h))
    # delta = (Xi ⊕ hip)\Xj # TODO -- THERE IS SOME MAJORLY SILLY ISSUE HERE
    # res[1:6] = veeEuler(delta)

  res[1:3] = (xi[1:3,idx] + [0,0,0.5]) - xj[1:3,idx] # cheating with 0.5 before rotation
  res[4:5] = xi[4:5,idx] - xj[4:5,idx]
  res[6] = wrapRad(wrapRad(xi[6,idx] + meas[1][idx]) - xj[6,idx])
  # sho1 = SE3( zeros(3), convert(SO3,so3([0,0,meas[1][idx]])) )
  # sho2 = SE3( [0,0,1.5], SO3(0) )
  #
  # del = (Xi ⊕ sho1 ⊕ sho2) \ Xj
  #
  # res[1:6] = veeEuler(del)
  nothing
end

mutable struct XJoint <: FunctorPairwise
  Zij::Distribution
end
function getSample(el::XJoint, N=1)
  return (rand(el.Zij, N),)
end
function (el::XJoint)(res, userdata, idx, meas, xi, xj)
  Xi = SE3(xi[1:3,idx],Euler(xi[4:6,idx]))
  Xj = SE3(xj[1:3,idx],Euler(xj[4:6,idx]))

  sho1 = SE3( zeros(3), convert(SO3,so3([meas[1][idx],0,0])) )
  sho2 = SE3( [0,0,1.0], SO3(0) )

  del = (Xi ⊕ sho1 ⊕ sho2) \ Xj

  res[1:6] = veeEuler(del)
  nothing
end

# setup visualization process and default drawings
vis = startdefaultvisualization()

N=100
fg = RoME.initfg()

# base
addVariable!(fg, :x1, dims=6)
pos = PriorPose3(MvNormal(zeros(6),1e-6*Matrix{Float64}(LinearAlgebra.I, 6,6)))
addFactor!(fg, [:x1], pos) # base
initVariable!(fg, :x1)

# torso
addVariable!(fg, :x2, dims=6)
hip = ZJoint(Normal(pi/3,0.1))
addFactor!(fg, [:x1, :x2], hip) # hio
initVariable!(fg, :x2)

# using IncrementalInference
#
# pts = evalFactor2(fg, fg.g.vertices[fg.fIDs[:x1x2]], fg.IDs[:x2])
#
# @show Base.mean(pts,2)

visualizeallposes!(vis, fg, drawtype=:max)

solveandvisualize(fg, vis)

Graphs.plot(fg.g)

visualizeDensityMesh!(vis, fg, :x2)

plotKDE(fg, :x2, dims=[4])

plotPose3Pairs(fg, :x2)

# torso
addVariable!(fg, :x3, dims=6)
should = XJoint(Normal(pi/4,0.1))
addFactor!(fg, [:x2, :x3], should) # hio
initVariable!(fg, :x3)

visualizeallposes!(vis, fg, drawtype=:max)

solveandvisualize(fg, vis)

# upper arm
addVariable!(fg, :x4, dims=6)
should = XJoint(Normal(pi/4,0.1))
addFactor!(fg, [:x3, :x4], should) # hio
initVariable!(fg, :x4)

visualizeallposes!(vis, fg, drawtype=:fit)

solveandvisualize(fg, vis)

plotPose3Pairs(fg, :x4)

# reachability Examples

N=300
fg = RoME.initfg()

# base
addVariable!(fg, :x1, dims=6)
pos = PriorPose3(MvNormal(zeros(6),1e-6*Matrix{Float64}(LinearAlgebra.I, 6,6)))
addFactor!(fg, [:x1], pos) # base
initVariable!(fg, :x1)

addVariable!(fg, :x2, dims=6)
hip = ZJoint(Uniform(-pi/3,pi/3))
addFactor!(fg, [:x1, :x2], hip) # hio
initVariable!(fg, :x2)

addVariable!(fg, :x3, dims=6)
should = XJoint(Uniform(-pi/4,pi/4))
addFactor!(fg, [:x2, :x3], should) # hio
initVariable!(fg, :x3)

visualizeallposes!(vis, fg, drawtype=:max)

solveandvisualize(fg, vis)

visualizeDensityMesh!(vis, fg, :x3)

plotPose3Pairs(fg, :x3)

getVal(fg, :x3)[5,1:10]

plotKDE(fg, :x3, dims=[6])

#

ext: jl


text:
#=
Search functions for Str strings

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on julia/base/strings/search.jl
=#

"""
    find(Fwd, pattern, string::AbstractString, start::Integer)

Find the next occurrence of `pattern` in `string` starting at position `start`.
`pattern` can be either a string, character, function, or regular expression
(in which case `string` must be of type `String`, or a `Str` type whose codeunit size is == 1,
or a `SubString` of one of those types)

The return value is a range of indexes where the matching sequence is found, such that
`s[find(Fwd, x, s, i)] == x`:

`find(Fwd, "substring", string, i)` = `start:end` such that
`string[start:end] == "substring"`, or `0:-1` if unmatched.

# Examples
```jldoctest
julia> find(Fwd, "z", "Hello to the world", 1)
0:-1

julia> find(Fwd, "o", "Hello to the world", 6)
8:8

julia> find(Fwd, "Julia", "JuliaLang", 2)
1:5

julia> find(Fwd, "z", "Hello to the world")
0:-1

julia> find(Fwd, "Julia", "JuliaLang")
1:5
```
"""
find(::Type{Fwd}, pat, str, pos)

"""
    find(First, pattern, string::AbstractString)

Find the first occurrence of `pattern` in `string` starting at position `start`.
`pattern` can be either a string, character, function, or regular expression
(in which case `string` must be of type `String`, or a `Str` type whose codeunit size is == 1,
or a `SubString` of one of those types)

The return value is a range of indexes where the matching sequence is found, such that
`s[find(Fwd, x, s, i)] == x`

`find(Fwd, "substring", string, i)` = `start:end` such that
`string[start:end] == "substring"`, or `0:-1` if unmatched.

# Examples
```jldoctest
julia> find(First, "z", "Hello to the world")
0:-1

julia> find(First, "z", "Hello to the world")
0:-1

julia> find(First, "Julia", "JuliaLang")
1:5
```
"""
find(::Type{First}, pat, str)

"""
    find(Rev, pattern::AbstractString, string::AbstractString, start::Integer)

Find the previous occurrence of `pattern` in `string` starting at position `start`.

The return value is a range of indexes where the matching sequence is found, such that
`s[find(Rev, x, s, i)] == x`

`find(Rev, "substring", string, i)` = `start:end` such that
`string[start:end] == "substring"`, or `0:-1` if unmatched.

# Examples
```jldoctest
julia> find(Rev, "z", "Hello to the world", 18)
0:-1

julia> find(Rev, "o", "Hello to the world", 18)
15:15

julia> find(Rev, "Julia", "JuliaLang", 6)
1:5
```
"""
find(::Type{Rev}, pat, str, pos)

"""
    find(Last, pattern::AbstractString, string::AbstractString)

Find the previous occurrence of `pattern` in `string`

The return value is a range of indexes where the matching sequence is found, such that
`s[find(First, x, s)] == x`

`find(First, "substring", string, i)` = `start:end` such that
`string[start:end] == "substring"`, or `0:-1` if unmatched.

# Examples
```jldoctest
julia> find(Last, "o", "Hello to the world")
15:15

julia> find(Last, "Julia", "JuliaLang")
1:5
```
"""
find(::Type{Last}, pat, str)

const _not_found = 0:-1

found(::Type{<:AbstractString}, v) = v != 0
find_result(::Type{<:AbstractString}, v) = v

nothing_sentinel(i) = first(i) == 0 ? nothing : i
Base.findfirst(a, b::Str)   = nothing_sentinel(find(First, a, b))
Base.findlast(a, b::Str)    = nothing_sentinel(find(Last, a, b))
Base.findnext(a, b::Str, i) = nothing_sentinel(find(Fwd, a, b, i))
Base.findprev(a, b::Str, i) = nothing_sentinel(find(Rev, a, b, i))
Base.findfirst(a::Str, b::AbstractString)   = nothing_sentinel(find(First, a, b))
Base.findlast(a::Str, b::AbstractString)    = nothing_sentinel(find(Last, a, b))
Base.findnext(a::Str, b::AbstractString, i) = nothing_sentinel(find(Fwd, a, b, i))
Base.findprev(a::Str, b::AbstractString, i) = nothing_sentinel(find(Rev, a, b, i))

function find(::Type{D}, fun::Function, str::AbstractString, pos::Integer) where {D<:Direction}
    pos < Int(D===Fwd) && (@boundscheck boundserr(str, pos); return 0)
    if pos > (len = ncodeunits(str))
        @boundscheck pos > len+1 && boundserr(str, pos)
        return 0
    end
    @inbounds is_valid(str, pos) || index_error(str, pos)
    _srch_pred(D(), fun, str, pos)
end

find(::Type{First}, fun::Function, str::AbstractString) =
    _srch_pred(Fwd(), fun, str, 1)
find(::Type{Last},  fun::Function, str::AbstractString) =
    _srch_pred(Rev(), fun, str, lastindex(str))

find(::Type{Next}, pat, str::AbstractString, pos) = find(Fwd, pat, str, pos)
find(::Type{Prev}, pat, str::AbstractString, pos) = find(Rev, pat, str, pos)

# Invert order of function and findop, to allow use of do blocks
find(pat::Function, ::Type{D}, str::AbstractString) where {D<:FindOp} = find(D, pat, str)
find(pat::Function, ::Type{D}, str::AbstractString, pos) where {D<:FindOp} = find(D, pat, str, pos)

find(::Type{D}, pred::P, str::AbstractString,
    pos::Integer) where {P<:Fix2{Union{typeof(==),typeof(isequal)}, <:AbsChar}, D<:Direction} =
    find(D, pred.x, str, pos)

function find(::Type{D}, ch::AbsChar, str::AbstractString, pos::Integer) where {D<:Direction}
    pos < Int(D===Fwd) && (@boundscheck boundserr(str, pos); return 0)
    if pos > (len = ncodeunits(str))
        @boundscheck pos > len+1 && boundserr(str, pos)
        return 0
    end
    # Only check if MultiCU
    is_multi(str) && (@inbounds is_valid(str, pos) || index_error(str, pos))
    # Check here if ch is valid for the type of string
    is_valid(eltype(str), ch) ? _srch_cp(D(), EncodingStyle(str), str, ch, pos, len) : 0
end

_get_dir(::Type{First}) = Fwd()
_get_dir(::Type{Last})  = Rev()

find(::Type{D}, ch::AbsChar, str::AbstractString) where {D<:Union{First,Last}} =
    ((len = ncodeunits(str)) == 0 || !is_valid(eltype(str), ch) ? 0
     : _srch_cp(_get_dir(D), EncodingStyle(str), str, ch,
                D === First ? 1 : lastindex(str), len))

function find(::Type{D}, needle::AbstractString, str::AbstractString,
              pos::Integer) where {D<:Direction}
    # Check for fast case of a single codeunit (should check for single character also)
    if (slen = ncodeunits(str)) == 0
        # Special case for empty string
        @boundscheck pos == 1 || pos == 0 || boundserr(str, pos)
        return ifelse(isempty(needle), 1:0, _not_found)
    end
    pos < 1 && (@boundscheck boundserr(str, pos) ; return _not_found)
    pos > slen && (@boundscheck pos > slen + 1 && boundserr(str, pos) ; return _not_found)
    @inbounds is_valid(str, pos) || index_error(str, pos)
    (tlen = ncodeunits(needle)) == 0 && return pos:pos-1
    (cmp = CanContain(str, needle)) === NoCompare() && return _not_found
    @inbounds ch, nxt = str_next(needle, 1)
    is_valid(eltype(str), ch) || return _not_found
    # Check if single character
    if nxt > tlen
        pos = _srch_cp(D(), EncodingStyle(str), str, ch, pos, slen)
        return pos == 0 ? _not_found : (pos:pos)
    end
    _srch_strings(D(), cmp, str, needle, ch, nxt, pos, slen, tlen)
end

function find(::Type{T}, needle::AbstractString, str::AbstractString) where {T<:Union{First,Last}}
    # Check for fast case of a single codeunit (should check for single character also)
    (slen = ncodeunits(str)) == 0 && return ifelse(isempty(needle), 1:0, _not_found)
    pos = T === First ? 1 : thisind(str, slen)
    (tlen = ncodeunits(needle)) == 0 && return pos:(pos-1)
    (cmp = CanContain(str, needle)) === NoCompare() && return _not_found
    @inbounds ch, nxt = str_next(needle, 1)
    is_valid(eltype(str), ch) || return _not_found
    # Check if single character
    if nxt > tlen
        pos = _srch_cp(_get_dir(T), EncodingStyle(str), str, ch, pos, slen)
        return pos:(pos - (pos == 0))
    end
    _srch_strings(_get_dir(T), cmp, str, needle, ch, nxt, pos, slen, tlen)
end

## Lower level search functions

@inline function _srch_pred(::Fwd, testf, str, pos)
    len = ncodeunits(str)
    while pos <= len
        @inbounds ch = str[pos]
        testf(ch) && return pos
        pos = nextind(str, pos)
    end
    0
end

@inline function _srch_pred(::Rev, testf, str, pos)
    while pos > 0
        @inbounds ch = str[pos]
        testf(ch) && break
        pos = prevind(str, pos)
    end
    pos
end

# AbstractString implementations of the generic find interfaces

@inline function _srch_codeunit(::Fwd, beg::Ptr{T}, cu::T, pos, len) where {T<:CodeUnitTypes}
    if sizeof(Cwchar_t) == sizeof(T) || T == UInt8
        pnt = _fwd_memchr(bytoff(beg, pos - 1), cu, len - pos + 1)
        pnt == C_NULL ? 0 : chrdiff(pnt, beg) + 1
    else
        beg -= sizeof(T)
        pnt = bytoff(beg, pos)
        fin = bytoff(beg, len)
        while pnt <= fin
            get_codeunit(pnt) == cu && return chrdiff(pnt, beg)
            pnt += sizeof(T)
        end
        0
    end
end

@inline function _srch_codeunit(::Rev, beg::Ptr{UInt8}, cu::UInt8, pos)
    pnt = _rev_memchr(beg, cu, pos)
    #println("beg=$beg, cu=$cu, pos=$pos => $pnt")
    pnt == C_NULL ? 0 : chrdiff(pnt, beg) + 1
end

@inline function _srch_codeunit(::Rev, beg::Ptr{T}, cu::T, pos) where {T<:Union{UInt16,UInt32}}
    pnt = bytoff(beg, pos)
    while (pnt -= sizeof(T)) >= beg && get_codeunit(pnt) != cu ; end
    #println("beg=$beg, cu=$cu, T=$T, pos=$pos => $pnt")
    chrdiff(pnt, beg) + 1
end

# _srch_cp is only called with values that are valid for that string type,
# and checking as already been done on the position (pos)
# These definitions only work for SingleCU types
_srch_cp(::Fwd, ::SingleCU, str::T, cp::AbsChar, pos, len) where {T<:Str} =
    @preserve str _srch_codeunit(Fwd(), pointer(str), cp%codeunit(T), pos, len)
_srch_cp(::Rev, ::SingleCU, str::T, cp::AbsChar, pos, len) where {T<:Str} =
    @preserve str _srch_codeunit(Rev(), pointer(str), cp%codeunit(T), pos)

function _srch_cp(::Fwd, cus, str, cp, pos, len)
    @inbounds while pos <= len
        str[pos] == cp && return pos
        pos = nextind(str, pos)
    end
    0
end

function _srch_cp(::Rev, cus, str, cp, pos, len)
    @inbounds while pos > 0
        str[pos] == cp && return pos
        pos = prevind(str, pos)
    end
    0
end

# Substring searching

# Todo: make use of CompareStyle trait to improve performance

"""Compare two strings, starting at nxtstr and nxtsub"""
@inline function _cmp_str(str, strpos, endpos, sub, subpos, endsub)
    while strpos <= endpos
        c, strnxt = str_next(str, strpos)
        d, subpos = str_next(sub, subpos)
        c == d || break
        subpos > endsub && return strpos
        strpos = strnxt
    end
    0
end

function _srch_strings(::Fwd, ::CompareStyle, str, needle, ch, subpos, pos, slen, tlen)
    cu = EncodingStyle(str)
    while (pos = _srch_cp(Fwd(), cu, str, ch, pos, slen)) != 0
        nxt = nextind(str, pos)
        res = _cmp_str(str, nxt, slen, needle, subpos, tlen)
        res == 0 || return pos:res
        pos = nxt
        pos > slen && break
    end
    _not_found
end

function _srch_strings(::Rev, ::CompareStyle, str, needle, ch, nxtsub, pos, slen, tlen)
    cu = EncodingStyle(str)
    while (prv = prevind(str, pos)) != 0 &&
        (loc = _srch_cp(Rev(), cu, str, ch, prv, slen)) != 0
        res = _cmp_str(str, nextind(str, loc), slen, needle, nxtsub, tlen)
        res == 0 || res > pos || return loc:res
        pos = loc
    end
    _not_found
end

_search_bloom_mask(ch) = UInt64(1) << (ch & 63)
_check_bloom_mask(mask, pnt, off) = (mask & _search_bloom_mask(get_codeunit(pnt, off))) == 0

function _srch_str_bloom(::Fwd, str, spnt, npnt, ch, pos, slen, tlen)
    tlast = get_codeunit(npnt, tlen)
    bloom_mask = _search_bloom_mask(tlast)
    skip = tlen - 1
    for j in 1:tlen-1
        cu = get_codeunit(npnt, j)
        bloom_mask |= _search_bloom_mask(cu)
        cu == tlast && (skip = tlen - j - 1)
    end

    pos -= 1
    while pos <= slen
        if get_codeunit(spnt, pos + tlen) == tlast
            pos += 1
            # check candidate
            j = 0
            while j < tlen - 1 && get_codeunit(spnt, pos + j) == get_codeunit(npnt, j + 1)
                j += 1
            end

            # match found
            j == tlen - 1 && return pos:_thisind(EncodingStyle(str), str, slen, spnt, pos + j)

            # no match, try to rule out the next character
            if pos <= slen && _check_bloom_mask(bloom_mask, spnt, pos + tlen)
                pos += tlen
            else
                pos += skip
            end
        elseif (pos += 1) <= slen && _check_bloom_mask(bloom_mask, spnt, pos + tlen)
            pos += tlen
        end
    end

    _not_found
end

function _srch_str_bloom(::Rev, str, spnt, npnt, ch, pos, slen, tlen)
    (slen < tlen || pos <= 0) && return 0

    tfirst = is_multi(str) ? get_codeunit(npnt) : codepoint(ch)
    bloom_mask = UInt64(0)
    skip = tlen - 1
    for j in skip:-1:0
        cu = get_codeunit(npnt, j)
        bloom_mask |= _search_bloom_mask(cu)
        cu == tfirst && j > 0 && (skip = j - 1)
    end

    pos = min(pos, slen) - tlen
    while pos >= 0
        if get_codeunit(spnt, pos) == tfirst
            # check candidate
            j = 1
            while j < tlen && get_codeunit(spnt, pos + j) == get_codeunit(npnt, j)
                j += 1
            end

            # match found
            j == tlen && return i

            # no match, try to rule out the next character
            pos -= ((pos > 0 && bloom_mask & _search_bloom_mask(get_codeunit(spnt, pos)) == 0)
                  ? tlen : skip)
        elseif pos > 0
            (bloom_mask & _search_bloom_mask(get_codeunit(spnt, pos)) == 0) && (pos -= tlen)
        end
        pos -= 1
    end

    0
end

# This should work for compatible CSEs, like ASCII & Latin, ASCII & UTF8, etc.
# See equals trait
_srch_strings(::Type{D}, ::Union{ByteCompare,WidenCompare}, str, needle,
              ch, nxtsub, pos, slen, tlen) where {D<:Direction} =
                  @preserve str needle _srch_str_bloom(D(), str, pointer(str), pointer(needle),
                                                       ch, pos, slen, tlen)

_occurs_in(needle, hay) = first(find(First, needle, hay)) != 0

# Avoid type piracy
occurs_in(needle::AbstractString, hay::Str) = _occurs_in(needle, hay)
occurs_in(needle::Str, hay::AbstractString) = _occurs_in(needle, hay)
occurs_in(needle::Chr, hay::AbstractString) = _occurs_in(needle, hay)
occurs_in(needle::Char, hay::Str)           = _occurs_in(needle, hay)
occurs_in(needle::Str, hay::Str)            = _occurs_in(needle, hay)

in(chr::Chr, str::AbstractString) = _occurs_in(chr, str)
in(chr::AbsChar,   str::Str)      = _occurs_in(chr, str)
in(pat::Str, str::AbstractString) = _occurs_in(pat, str)
in(pat::AbstractString, str::Str) = _occurs_in(pat, str)
in(pat::Str, str::Str)            = _occurs_in(pat, str)


text:
# emulating exmaple found on pardiso_unsym_schur.cpp found on https://pardiso-project.org

n = 16

ia = [0, 9, 17, 24, 31, 38, 46, 53, 61, 63, 65, 67, 69, 71, 73, 75, 77] .+ 1

ja = [0, 2, 5, 6, 8, 9, 10, 11, 12, 1, 2, 4, 8, 9, 10, 11, 12, 2, 7, 8, 9,
                     10, 11, 12, 3, 6, 8, 9, 10, 11, 12, 1, 4, 8, 9, 10, 11, 12, 2, 5, 7,
                     8, 9, 10, 11, 12, 1, 6, 8, 9, 10, 11, 12, 2, 6, 7, 8, 9, 10, 11, 12,
                     0, 8, 1, 9, 2, 10, 3, 11, 4, 12, 5, 13, 6, 14, 7, 15] .+ 1

a = [
    complex(7.0000000000000000, 1.),
    complex(1.0000000000000000, 1.),
    complex(2.0000000000000000, 1.),
    complex(7.0000000000000000, 1.),
    complex(0.8813528146570788, 0.),
    complex(0.9693082988664113, 0.),
    complex(0.0063764102307857, 0.),
    complex(0.0656724982252748, 0.),
    complex(0.1285133482338328, 0.),
    complex(-4.0000000000000000, 0.),
    complex(8.0000000000000000, 1.),
    complex(2.0000000000000000, 1.),
    complex(0.0246814022448006, 0.),
    complex(0.3913069588115609, 0.),
    complex(0.9199013631348676, 0.),
    complex(0.2300406237580992, 0.),
    complex(0.6371356638478847, 0.),
    complex(1.0000000000000000, 1.),
    complex(5.0000000000000000, 1.),
    complex(0.3412064600197927, 0.),
    complex(0.3135851048672553, 0.),
    complex(0.0179622413055945, 0.),
    complex(0.1170403198125443, 0.),
    complex(0.7464881699547887, 0.),
    complex(7.0000000000000000, 0.),
    complex(9.0000000000000000, 1.),
    complex(0.4214056775317185, 0.),
    complex(0.5533135799374451, 0.),
    complex(0.0294164798770650, 0.),
    complex(0.9897710614201980, 0.),
    complex(0.8052672976354331, 0.),
    complex(-4.0000000000000000, 1.),
    complex(0.0000000000000000, 0.),
    complex(0.0800111586678662, 0.),
    complex(0.7920480624973907, 0.),
    complex(0.7114519163682994, 0.),
    complex(0.0705164910828937, 0.),
    complex(0.2690271315860681, 0.),
    complex(7.0000000000000000, 1.),
    complex(3.0000000000000000, 1.),
    complex(8.0000000000000000, 0.),
    complex(0.0793538091733186, 0.),
    complex(0.7983140830802098, 0.),
    complex(0.5384141865401406, 0.),
    complex(0.0600236559152354, 0.),
    complex(0.4339797789833593, 0.),
    complex(1.0000000000000000, 1.),
    complex(11.0000000000000000, 1.),
    complex(0.0723762491986621, 0.),
    complex(0.8632585405599013, 0.),
    complex(0.5245885805189759, 0.),
    complex(0.6615771015291297, 0.),
    complex( 0.4018064872641355, 0.),
    complex(-3.0000000000000000, 1.),
    complex( 2.0000000000000000, 1.),
    complex( 5.0000000000000000, 0.),
    complex( 0.9003410060810889, 0.),
    complex( 0.7979622044473094, 0.),
    complex( 0.5022399304161738, 0.),
    complex( 0.3440577694777694, 0.),
    complex(0.6081423394349520, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.),
    complex(1.0000000000000000, 0.),
    complex(0.0000000000000000, 0.)
    ]

M = SparseMatrixCSC(n,n,ia,ja,a)

ps = PardisoSolver()
set_matrixtype!(ps,Pardiso.COMPLEX_NONSYM)
S = schur_complement(ps,M,8)

m = real(M)
set_matrixtype!(ps,Pardiso.REAL_NONSYM)
s = schur_complement(ps,m,8)

A = M[1:n-8,1:n-8]
A⁻¹ = inv(Matrix(A))
rA⁻¹ = inv(real(Matrix(A)))
B = M[1:n-8,n-7:end]
C = M[n-7:end,1:n-8]
D = M[n-7:end,n-7:end]


text:
using CanonicalTraits
using Test
import LinearAlgebra

@trait Monoid{A} begin
    mempty :: Type{A} => A
    (⊕)    :: [A, A] => A
end

@implement Monoid{Int} begin
    mempty(_) = 0
    ⊕         = +
end

@implement Monoid{Vector{T}} where T begin
    mempty(_) = T[]
    ⊕         = vcat
end

@trait ToString{A} begin # like Show in Haskell
    to_str :: A => String
    to_str = Base.string # default implementation
end

@implement ToString{A} where A # for all types

# 100 ⊕ 2 |> println
# @btime 100 ⊕ 2
# @btime 100 + 2
# @btime vcat([1, 2, 3], [3, 4, 5])
# @btime [1, 2, 3] ⊕ [3, 4, 5]
# @btime 1 ⊕ 2

@testset "simple" begin
    # Write your own tests here.
    @test_throws Any begin
        "2" ⊕ "3"
    end

    @test 102 == begin
        100 ⊕ 2
    end

    @test [1, 2, 3, 3] == begin
        [1, 2] ⊕ [3, 3]
    end

    @test map(to_str, ["123", 10]) == map(string, ["123", 10])
end

using Polynomials

function vect_infer_helper
end

@trait Vect{F, V} where {F = vect_infer_helper(V)} begin
    scalar_mul :: [F, V] => V
    scalar_div :: [V, F] => V

    vec_add    :: [V, V] => V
    vec_sub    :: [V, V] => V

    scalar_add :: [F, V] => V
    scalar_sub :: [V, F] => V

    scalar_div(vec :: V, scalar :: F) = scalar_mul(one(F)/scalar, vec)
    scalar_sub(vec :: V, scalar :: F) = scalar_add(-scalar, vec)
    vec_sub(vec1 :: V, vec2 :: V)     = vec_add(vec1, scalar_mul(-one(F), vec2))
end

vect_infer_helper(::Type{Polynomial{T}}) where T = T

@implement Vect{F, Polynomial{F}} where F <: Number begin
    function scalar_mul(num::F, vec::Polynomial{F}) where F <: Number
        num * vec
    end
    function vec_add(vec1 :: Polynomial{F}, vec2 :: Polynomial{F}) where F <: Number
        vec1 + vec2
    end
    function scalar_add(scalar::F, vec::Polynomial{F}) where F <: Number
        scalar + vec
    end
end

@trait Vect{F, V} >: Dot{F <: Number, V} where {F = vect_infer_helper(V)} begin
    dot :: [V, V] => F
    gram_schmidt :: [V, Vector{V}] => V
    function gram_schmidt(v :: V, vs :: Vector{V})::V where F <: Number
        for other in vs
            coef = dot(v, other) / dot(other, other)
            v = vec_sub(v, scalar_mul(coef, other))
        end
        scalar_div(v, sqrt(dot(v, v)))
    end
end

@implement Dot{F, Polynomial{F}} where F <: Number begin
    function dot(v1 :: Polynomial{F}, v2 :: Polynomial{F})::Real where F <: Number
            f = Polynomials.integrate(v1 * v2)
            f(1) - f(-1)
    end
end

vect_infer_helper(::Type{Tuple{F, F}}) where F<:Number = F

@implement Vect{F, Tuple{F, F}} where F <: Number begin
    scalar_add(num, vec) =
        (vec[1] + num, vec[2] + num)
    vec_add(vec1, vec2) =
        (vec1[1] + vec2[1], vec1[2] + vec2[2])
    scalar_mul(num, vec) =
        (num * vec[1], num * vec[2])
end

@implement Dot{F, Tuple{F, F}} where F <: Number begin
    function dot(v1, v2)
        LinearAlgebra.dot(F[v1[1], v1[2]], F[v2[1], v2[2]])
    end
end

@testset "polynomial orthogonalization" begin

    @test scalar_add(5.0, Polynomial([2.0, 1.0])) == Polynomial([7.0, 1.0])
    fx1 = Polynomial([1.0])
    fx2 = Polynomial([0.0, 1.0])
    T = typeof(fx1)
    fx1_ot = gram_schmidt(fx1, T[])
    fx2_ot = gram_schmidt(fx2, T[fx1_ot])

    @test dot(fx1_ot, fx2_ot) ≈ 0
    @test dot(fx1_ot, fx1_ot) ≈ 1
    @test dot(fx2_ot, fx2_ot) ≈ 1

    fx1 = (1.0, 2.0)
    fx2 = (3.0, 5.0)
    T = typeof(fx1)
    fx1_ot = gram_schmidt(fx1, T[])
    fx2_ot = gram_schmidt(fx2, T[fx1_ot])

    @test dot(fx1_ot, fx2_ot) + 1.0 ≈ 1.0
    @test dot(fx1_ot, fx1_ot) ≈ 1.0
    @test dot(fx2_ot, fx2_ot) ≈ 1.0

end

using MLStyle
function type_constructor_from_hkt end
function type_argument_from_hkt end
function type_app end

# type app representation
struct App{Cons, K₀}
    injected :: Any
end

@trait Higher{Cons, K₀, K₁} where {
    Cons=type_constructor_from_hkt(K₁),
    K₀=type_argument_from_hkt(K₁),
    K₁=type_app(Cons, K₀)
} begin
    inj :: K₁ => App{Cons, K₀}
    inj(data::K₁) = App{Cons, K₀}(data)
    prj :: App{Cons, K₀} => K₁
    prj(data::App{Cons, K₀})::K₁ = data.injected
end

abstract type HKVect end
Base.@pure type_constructor_from_hkt(::Type{Vector{T}}) where T = HKVect
Base.@pure type_argument_from_hkt(::Type{Vector{T}}) where T = T
Base.@pure type_app(::Type{HKVect}, ::Type{T}) where T = Vector{T}
@implement Higher{HKVect, T, Vector{T}} where T

@testset "higher kinded" begin
    hkt_vect = inj([1, 2, 3])
    @test (hkt_vect |> typeof) == App{HKVect, Int}
    @test prj(hkt_vect) == [1, 2, 3]
end

@trait P{A} begin
    fx :: A => Int
end

@implement P{Symbol} begin
    fx(x) = 1
end

@implement P{Tuple{T, T}} where T begin
    fx(x) = fx(x[1]) + fx(x[2])
end

@testset "mutually referencing" begin
    @test fx(:a) == 1
    @test fx((:a, :b)) == 2
end

@trait Add1{T <: Number} begin
    add1 :: [T] => T
end

@trait Add1{T} >: Addn{T <: Number} begin
    addn :: [Int, T] => T
    addn(n, x) = let s = x; for i in 1:n; s = add1(s) end; s; end
end

@implement Add1{Int} begin
    add1(x) = x + 1
end

@implement Addn{Int}

@testset "class inheritance" begin
    @test add1(1) == 2
    @test addn(5, 1) == 6
    @test "Not implemented trait Add1 for (Float64)." == try
        addn(2, 1.9)
        ""
    catch e
        strip(e.msg)
    end
end

@implement! Add1{T} >: Add1{Vector{T}} where T begin
    add1(xs) = add1.(xs)
end

@testset "instance inheritance" begin
    @test add1([1, 2, 3]) == [2, 3, 4]
    @test "Not implemented trait Add1 for (Float64)." == try
        add1([1.])
    catch e
        strip(e.msg)
    end
end


text:
"""
    groupby(f, xs)
group values `x ∈ xs` by the result of applying `f`.
Returns a `Dict` with keys `yi` and values `[xi1, xi2,...]` such that `f(xij) = yi`.

# Example
```julia-repl
julia> groupby(isodd,1:10)
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [2, 4, 6, 8, 10]
  true  => [3, 5, 7, 9]
```
"""
function groupby(f::Function, xs)
    # maintain order in xs
    sxs = Iterators.Stateful(xs)
    T = eltype(sxs)
    x = first(sxs)
    y = f(x)
    dict = Dict{typeof(y), Vector{T}}()
    dict[y] = T[x]
    for x in sxs
        y = f(x)
        haskey(dict, y) || (dict[y] = T[])
        push!(dict[y], x)
    end
    return dict
end

"""
    gatherby(f, xs)
Like `groupby` but only returns values, i.e. elements of `xs` in `Vector`s such that
`f` applied to an element of a group is the same as for any other element of that
`Vector`.

# Example
```julia-repl
julia> gatherby(isodd, 1:10)
2-element Array{Array{Int64,1},1}:
 [2, 4, 6, 8, 10]
 [3, 5, 7, 9]
```
"""
gatherby(f::Function, xs) = collect(values(groupby(f,xs)))


text:
using VersionParsing
using Test

@test vparse("3.7.2a4") == v"3.7.2-a4" == vparse("version 3.7.2a4: the best version")
@test vparse("2.1.0-python3_5") == v"2.1.0-python35" # Plots.jl#1432
@test vparse("0.99.1.1") == v"0.99.1+1" # julia#7282
@test vparse("2.1.0.post806+g905465b") == v"2.1.0+post806.g905465b" # from julia#7282

# examples from npm docs (https://docs.npmjs.com/misc/semver)
@test_throws ArgumentError vparse("a.b.c")
@test vparse("  =v1.2.3   ") == v"1.2.3"
@test vparse("v2") == vparse("version 2.0") == v"2.0.0"
@test vparse("42.6.7.9.3-alpha") == v"42.6.7-alpha+9.3"

@test vparse("1.2.3.8-4.2") == v"1.2.3-4.2+8"
@test vparse("1.2.3.8-4..b2+9.7++2-0") == v"1.2.3-4.b2+8.9.7.2-0"

# Debian-style epoch:upstream version numbers
@test vparse("4:4.7.4-0ubuntu8") == vparse("0:4.7.4-0ubuntu8") == v"4.7.4-0ubuntu8"

@test vparse(".3") == v"0.3.0"
@test vparse("2.1.0...+3") == v"2.1.0+3"

@test vparse("3,1,2") == v"3.1.2" == vparse("version 3.1.2, a good version")

@test vparse("A darn good version. 1.2") == v"1.2.0"
@test vparse("A darn good version. .2") == v"0.2.0"
@test_throws ArgumentError vparse("..2")


text:
module HTMLSanitizer

using Gumbo

export sanitize

"""
    sanitize(input::AbstractString; isfragment = true, whitelist = WHITELIST, prettyprint = false)

Sanitizes the HTML input according to `whitelist`.

- `isfragment`: If true, removes enclosing `<HTML>` tags from the output.
- `whitelist`: Whitelist for allowed elements and attributes.
- `prettyprint`: Returns a prettier multiline string instead of a somewhat minified version.
"""
function sanitize(input::AbstractString; isfragment = true, whitelist = WHITELIST, prettyprint = false)
    doc = parsehtml(input, preserve_whitespace=true)
    sanitize_bfs(doc.root, whitelist)

    out = IOBuffer()
    print(out, doc.root, pretty = prettyprint)

    out = String(take!(out))

    if isfragment
        out = replace(out, r"^<HTML>" => "")
        out = replace(out, r"</HTML>$" => "")
    else
        return out
    end
end

reparent!(node, parent) = node.parent = parent

function sanitize_bfs(tree, whitelist)
    i = 1
    while i <= length(tree.children)
        el = tree.children[i]

        sanitized = sanitize_element(el, whitelist)
        if sanitized isa Vector
            # reparent all nodes
            reparent!.(sanitized, Ref(tree))
            splice!(tree.children, i, sanitized)
            # don't increment i here so the newly inserted nodes are sanitized in the next iteration
        else
            # reparent node
            reparent!(sanitized, tree)
            tree.children[i] = sanitized
            i += 1
        end
    end
    sanitize_bfs.(tree.children, Ref(whitelist))
end

sanitize_bfs(tree::HTMLText, whitelist) = nothing

function sanitize_element(el::HTMLElement{TAG}, whitelist) where TAG
    tag = string(TAG)

    @debug("Sanitizing `$(tag)`.")

    if !(tag in get(whitelist, :elements, []))
        @debug("Element `$(tag)` not in whitelist.")
        if tag in get(whitelist, :remove_contents, [])
            @debug("Removing contents for `$(tag)`.")
            return Gumbo.HTMLText("")
        end
        @debug("Replacing `$(tag)` with its contents.")
        out = el.children
        return isempty(out) ? Gumbo.HTMLText("") : out
    end

    el = sanitize_attributes(el, whitelist)

    return el
end

sanitize_element(el::HTMLElement{:HTML}, whitelist) = el

sanitize_element(el::HTMLText, whitelist) = el

const REGEX_PROTOCOL = r"\A\s*([^\/#]*?)(?:\:|&#0*58|&#x0*3a)"i

sanitize_attributes(el, whitelist) = el

function sanitize_attributes(el::HTMLElement{TAG}, whitelist) where TAG
    tag = string(TAG)
    attributes = attrs(el)
    protocols = get(get(whitelist, :protocols, Dict()), tag, Dict())

    attributes_for_tag = get(get(whitelist, :attributes, Dict()), tag, [])
    attributes_for_all = get(get(whitelist, :attributes, Dict()), :ALL, [])

    for (attr, val) in attributes
        if !(attr in attributes_for_tag) && !(attr in attributes_for_all)
            # not in whitelist, so remove the attribute altogether
            @debug("Deleting attribute `$(attr)` in element `$(tag)` (not in whitelist).")
            delete!(attributes, attr)
        elseif haskey(protocols, attr)
            # allowed, but only specific values are ok
            is_acceptable = false

            if occursin(REGEX_PROTOCOL, val)
                # looks like a protocol is specified
                if any(startswith.(Ref(lowercase(val)), string.(protocols[attr])))
                    is_acceptable = true
                end
            else
                if :relative in protocols[attr] && is_relative_url(val)
                    is_acceptable = true
                end
            end

            if !is_acceptable
                @debug("Deleting attribute `$(attr)` in element `$(tag)` (does not conform to protocol).")
                delete!(attributes, attr)
            end
        end
    end

    return el
end

# A relative URL either
# 1. starts with `/` (root-relative).
# 2. starts with `//` (protocol-relative).
# 3. starts with `../`/`./` (relative directory traversal)
# 4. doesn't start with either of the above and doesn't start with a protocol (e.g. `foo/bar.html`)
function is_relative_url(url)
    if occursin(r"^\.?\.?//?"i, url)
        return true
    else
        return !occursin(r"^\w+://"i, url)
    end
end

"""
Default whitelist. Allows many elements and attributes, but crucially removes `<script>` elements
as well as `style` attributes.
"""
const WHITELIST = Dict(
    :elements => [
        "h1","h2","h3","h4","h5","h6","h7","h8","br","b","i","strong","em","a","pre","code","img","tt",
        "div","ins","del","sup","sub","p","ol","ul","table","thead","tbody","tfoot","blockquote",
        "dl","dt","dd","kbd","q","samp","var","hr","ruby","rt","rp","li","tr","td","th","s","strike",
        "summary","details","caption","figure","figcaption","abbr","bdo","cite","dfn","mark",
        "small","span","time","wbr","center"
    ],
    :remove_contents => ["script"],
    :attributes => Dict(
        "a"          => ["href"],
        "img"        => ["src", "longdesc"],
        "div"        => ["itemscope", "itemtype"],
        "blockquote" => ["cite"],
        "del"        => ["cite"],
        "ins"        => ["cite"],
        "q"          => ["cite"],
        :ALL         => [
            "abbr", "accept", "accept-charset",
            "accesskey", "action", "align", "alt",
            "aria-describedby", "aria-hidden", "aria-label", "aria-labelledby",
            "axis", "border", "cellpadding", "cellspacing", "char",
            "charoff", "charset", "checked",
            "clear", "cols", "colspan", "color",
            "compact", "coords", "datetime", "dir",
            "disabled", "enctype", "for", "frame",
            "headers", "height", "hreflang",
            "hspace", "ismap", "label", "lang",
            "maxlength", "media", "method",
            "multiple", "name", "nohref", "noshade",
            "nowrap", "open", "prompt", "readonly", "rel", "rev",
            "rows", "rowspan", "rules", "scope",
            "selected", "shape", "size", "span",
            "start", "summary", "tabindex", "target",
            "title", "type", "usemap", "valign", "value",
            "vspace", "width", "itemprop"
        ]
    ),
    :protocols => Dict(
        "a"          => Dict("href" => ["http", "https", "mailto", :relative]),
        "img"        => Dict(
            "src"      => ["http", "https", :relative],
            "longdesc" => ["http", "https", :relative]),
        "blockquote" => Dict("cite" => ["http", "https", :relative]),
        "del"        => Dict("cite" => ["http", "https", :relative]),
        "ins"        => Dict("cite" => ["http", "https", :relative]),
        "q"          => Dict("cite" => ["http", "https", :relative]),
    )
)

"""
Similar to the default whitelist, but only allows very few elements types.
"""
const LIMITED = merge(WHITELIST, Dict(
    :elements => ["b", "i", "strong", "em", "a", "pre", "code", "img", "ins", "del", "sup", "sub", "mark", "abbr", "p", "ol", "ul", "li"]
))

end # module


text:
using BenchmarkTools
println("\n# levendist(kitten, sitting)")
s, t = "kitten", "sitting"
println(" - Recursive:")
@btime levendist(s, t)
println(" - Iterative:")
@btime levendist1(s, t)
println("\n# levendist(rosettacode, raisethysword)")
s, t = "rosettacode", "raisethysword"
println(" - Recursive:")
@btime levendist(s, t)
println(" - Iterative:")
@btime levendist1(s, t)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct ZygoteRuleConfig{CTX<:AContext} <: RuleConfig{Union{HasReverseMode,NoForwardsMode}}
  context::CTX
end
ZygoteRuleConfig() = ZygoteRuleConfig(Context())

_is_rrule_redispatcher(m::Method) = m.sig == Tuple{typeof(rrule), RuleConfig, Vararg}

"""
  has_chain_rrule(T)

For a type-tuple `T` e.g. `Tuple{typeof(f), Int, Float64}`, checks if there is a `rrule` defined for it.
Excluding the generic fallback.
The first return value is `true` if the `rrule` exists, `false` otherwise.
If it does not, then the second argument is a list of edges to attach to the CodeInfo for a generated function,
such that if a suitable rule is defined later, the generated function will recompile.
"""
function has_chain_rrule(T)
  config_T, arg_Ts = Iterators.peel(T.parameters)
  configured_rrule_m = meta(Tuple{typeof(rrule), config_T, arg_Ts...})
  if _is_rrule_redispatcher(configured_rrule_m.method)
    # The config is not being used:
    # it is being redispatched without config, so we need the method it redispatches to
    rrule_m = meta(Tuple{typeof(rrule), arg_Ts...})
    # Thus any no_rrule that might apply must also not have a config because if there was a
    # no_rrule with a config that applied then there would also be a rrule with config that applied
    no_rrule_m = meta(Tuple{typeof(ChainRulesCore.no_rrule), arg_Ts...})
  else
    # Not being redispatched: it does have a config
    rrule_m = configured_rrule_m
    # Thus any no_rrule that might apply must also have a config because if it applied
    # it will be identical, and if it doesn't we don't care what it is.
    no_rrule_m = meta(Tuple{typeof(ChainRulesCore.no_rrule), config_T, arg_Ts...})
  end

  # To understand why we only need to check if the sigs match between no_rrule_m and rrule_m
  # in order to decide if to use, one must consider the following facts:
  # - for every method in `no_rrule` there is a identical one in `rrule` that returns nothing
  # - this includes the general fallback `rrule(::Any...)=nothing`.
  # - a configured rrule/no_rrule is always more specific than a otherwise equivalent unconfigured rrule/no_rrule
  #
  # Consider the following truth table, for what can occur:
  # rrule: fallback, no_rrule: fallback =>  matches => do not use rrule.
  # rrule: specific, no_rrule: fallback => !matches => do use rrule, as haven't opted out.
  # rrule: fallback, no_rrule: specific =>  IMPOSSIBLE, every no_rule is identical to some rrule
  # rrule: specific, no_rrule: specific =>  matches => do not use rrule as opted out
  # rrule: specific, no_rrule: general  => !matches => do use rrule as a more specific rrule takes preciedent over more general opted out
  # rrule: general , no_rrule: specific =>  IMPOSSIBLE, every no_rule us identical to some rrule so can't have a more general rrule being hit, as the specific one would hit first
  #
  # Note that the fallback cases are the same outcome as the general cases as fallback is just most general.
  # It can be seen that checking if it matches is the correct way to decide if we should use the rrule or not.

  do_not_use_rrule = matching_cr_sig(no_rrule_m, rrule_m)
  if do_not_use_rrule
    # Return instance for configured_rrule_m as that will be invalidated
    # directly if configured rule added, or indirectly if unconfigured rule added
    # Do not need an edge for `no_rrule` as no addition of methods to that can cause this
    # decision to need to be revisited (only changes to `rrule`), since we are already not
    # using the rrule, so not using more rules wouldn't change anything.
    return false, configured_rrule_m.instance
  else
    # Otherwise found a rrule, no need to add any edges for `rrule`, as it will generate
    # code with natural edges if a new method is defined there.
    # We also do not need an edge to `no_rrule`, as any time a method is added to `no_rrule`
    # a corresponding method is added to `rrule` (to return `nothing`), thus we will already
    # be revisiting this decision when a new opt-out is added.
    return true, nothing
  end
end

matching_cr_sig(t, s) = matching_cr_sig(t.method.sig, s.method.sig)
matching_cr_sig(::DataType, ::UnionAll) = false
matching_cr_sig(::UnionAll, ::DataType) = false
matching_cr_sig(t::Type, s::Type) = type_tuple_tail(t) == type_tuple_tail(s)

type_tuple_tail(d::DataType) = Tuple{d.parameters[2:end]...}
function type_tuple_tail(d::UnionAll)
    body = Base.unwrap_unionall(d)
    body_tt = type_tuple_tail(body)
    return Base.rewrap_unionall(body_tt, d)
end

"""
    is_kwfunc(sigt...)

Determines if `sigt` is the type signature of a kwfunction.
Each element of `sigt` should be a type.
Either the first 3 types are a kwfunc type, a NamedTuple and the matching base function type,
or the first argument is the base function type and it is not a kwfunction.
the remaining types in `sigt` are the types of the argument.

"""
is_kwfunc(::Vararg{Any}) = false
is_kwfunc(k, ::Type{<:NamedTuple}, f, args...) = k===Core.kwftype(f)

"""
    wrap_chainrules_output(x)

Convert `x` from the differentials types ChainRules uses to the format Zygote uses internally.
"""
@inline wrap_chainrules_output(x) = x
@inline wrap_chainrules_output(x::AbstractThunk) = wrap_chainrules_output(unthunk(x))  # For now we are just not going to deal with thunks
@inline wrap_chainrules_output(x::Tuple) = map(wrap_chainrules_output, x)
# Zygote convention: even if many AbstractZero partials (i.e. multi-input function), make just 1 nothing.
@inline wrap_chainrules_output(x::Tuple{Vararg{ChainRules.AbstractZero}}) = nothing
@inline wrap_chainrules_output(x::ChainRules.AbstractZero) = nothing
for T_outer in (:Tuple, :NamedTuple)
  # we create separate methods rather than using a `Union` + an `if` so that we avoid a
  # branch that changes output type, because nested AD on that kinda thing makes Zygote less
  # than happy.
  @eval @inline function wrap_chainrules_output(x::ChainRules.Tangent{P, T}) where {P, T<:$T_outer}
    xp = map(wrap_chainrules_output, canonicalize(x))
    ChainRulesCore.backing(xp)  # this is accessing ChainRulesCore internals, but it is prob safe enough, and it is fastest
  end
end
wrap_chainrules_output(dxs::AbstractArray{<:Number}) = dxs
wrap_chainrules_output(dxs::AbstractArray{<:AbstractArray{<:Number}}) = dxs
wrap_chainrules_output(dxs::AbstractArray) = map(wrap_chainrules_output, dxs)
#=
# As an optimisation, we can convert by `reinterpret` for bitstypes, e.g. arrays of tuples of numbers
@inline function wrap_chainrules_output(dxs::AbstractArray{<:ChainRules.Tangent{<:Any, B}}) where {B}
  if isbitstype(B)
    # B is the backing type. It still contains NoTangent etc, which need converting to Nothing
    reinterpret(wrap_chainrules_output(B), dxs)
  else
    map(wrap_chainrules_output, dxs)
  end
end
wrap_chainrules_output(::Type{<:AbstractZero}) = Nothing
wrap_chainrules_output(::Type{NamedTuple{L,T}}) where {L,T} = NamedTuple{L,wrap_chainrules_output(T)}
@generated function wrap_chainrules_output(::Type{T}) where T<:Tuple
  inner = map(wrap_chainrules_output, T.parameters)
  :(Tuple{$(inner...)})
end
=#

"""
    wrap_chainrules_input(dx)

Convert `dx` from the format Zygote uses internally to differentials types ChainRules uses.
"""
@inline wrap_chainrules_input(dx) = dx
@inline wrap_chainrules_input(::Nothing) = ChainRules.ZeroTangent()
@inline wrap_chainrules_input(::Tuple{Vararg{Nothing}}) = ChainRules.ZeroTangent()
@inline wrap_chainrules_input(::AbstractArray{Nothing}) = ChainRules.ZeroTangent()
@inline function wrap_chainrules_input(dxs::Union{Tuple, NamedTuple})
  xp = map(wrap_chainrules_input, dxs)
  # This produces Tangent{Any} since it does not get to see the primal, `x`.
  ChainRulesCore.Tangent{Any, typeof(xp)}(xp)
end
# For mutable types, including x=Ref(1), Zygote makes Ref{Any}(::NamedTuple)
@inline wrap_chainrules_input(dx::Ref) = wrap_chainrules_input(dx[])
# For arrays, whitelist the safe ones, but always look inside Any[]:
@inline wrap_chainrules_input(dxs::AbstractArray{<:Number}) = dxs
@inline wrap_chainrules_input(dxs::AbstractArray{<:AbstractArray{<:Number}}) = dxs
@inline wrap_chainrules_input(dxs::AbstractArray) = map(wrap_chainrules_input, dxs)

#=
# Could `reinterpret` instead here? See issue 1112.
# One easy case, might be this:
@inline wrap_chainrules_input(xs::Base.ReinterpretArray{<:NamedTuple, <:Tangent}) = parent(xs)

# This is for `z2d` reinterpret below:
wrap_chainrules_input(::Type{Nothing}) = NoTangent
wrap_chainrules_input(::Type{NamedTuple{L,T}}) where {L,T} = NamedTuple{L,wrap_chainrules_input(T)}
@generated function wrap_chainrules_input(::Type{T}) where T<:Tuple
  inner = map(wrap_chainrules_input, T.parameters)
  :(Tuple{$(inner...)})
end
=#

"""
  _project(x, dx)

Uses `ChainRulesCore.ProjectTo` to standardise the gradient `dx` for type & shape.
Also handles some Zygote-specific corrections, such as `x::Array, dx::Tuple`.
Safe to apply to arbitrary input.
"""
@inline function _project(x, dx)
  wrap_chainrules_output(ProjectTo(x)(zygote2differential(dx, x)))
end

# Restore splatted arrays
_project(x::AbstractArray, dx::Tuple) = _project(x, reshape(collect(dx), axes(x)))

# Piracy:
# CRC likes Tangent{AbstractArray}, but Zygote makes Tangent{Any}
# in particular this would hit https://github.com/JuliaDiff/ChainRulesCore.jl/blob/2ec2549b73b22bc08f554dae864fb650cfb9c3d7/src/projection.jl#L139
# if we were not losing track of the Primal in the Tangent
# This type piracy is just giving up that safety check.
(project::ProjectTo{AbstractArray})(dx::Tangent) = dx

"""
  ZBack{F}(back) <: Function

Wrapper for a ChainRules pullback `back`, that causes it to follow Zygote conventions.
(A functor here is used rather than a closure to avoid boxing issues);
"""
struct ZBack{F} <: Function
  back::F
end
@inline (s::ZBack)(dy) = wrap_chainrules_output(s.back(wrap_chainrules_input(dy)))
# `nothing->nothing` can be deleted after https://github.com/FluxML/Zygote.jl/issues/603
# though it might be worth keeping as a performance optimization (benchmarking pending)
@inline (s::ZBack)(::Nothing) = nothing

"""
    chain_rrule(config, f, args...)

Returns a the (primal) value of `f(args...)` and a pullback, by invoking `ChainRulesCore.rrule(f, args...)`.
The pullback is appropriately wrapped up to follow Zygote conventions.
"""
@inline function chain_rrule(config, f, args...)
  y, back = rrule(config, f, args...)
  return y, ZBack(back)
end

"""
  chain_rrule_kw(config, kwf, kwargs, f, args...)

As per [`chain_rrule`](@ref) but with support for kwargs.
`kwf` should be the kwfunc matching to `f`, and `kwargs` are a `NamedTuple` of keyword arguments.
"""
@inline function chain_rrule_kw(config, kwf, kwargs, f, args...)
  y, back = rrule(config, f, args...; kwargs...)
  function kw_zpullback(dy)
    dxs = ZBack(back)(dy)
    if dxs === nothing  # if dxs is nothing, then all partiaols are nothing
      # Zygote convention is a single nothing no mather how partials, if all are nothing
      return nothing
    else
      return (nothing, nothing, dxs...)  # first two nothings are for kwfunc and kwargs
    end
  end
  return y, kw_zpullback
end

function ChainRulesCore.rrule_via_ad(config::ZygoteRuleConfig, f_args...; kwargs...)
    # first check whether there is an `rrule` which handles this directly
    direcct = rrule(config, f_args...; kwargs...)
    direcct === nothing || return direcct

    # create a closure to work around _pullback not accepting kwargs
    # but avoid creating a closure unnecessarily (pullbacks of closures do not infer)
    y, pb = if !isempty(kwargs)
        kwf() = first(f_args)(Base.tail(f_args)...; kwargs...)
        _y, _pb = _pullback(config.context, kwf)
        _y, Δ -> first(_pb(Δ)).f_args  # `first` should be `only`
    else
        _pullback(config.context, f_args...)
    end

    ad_pullback(Δ) = zygote2differential(pb(wrap_chainrules_output(Δ)), f_args)
    return y, ad_pullback
end

"""
    zygote2differential(dx, primal)

Convert input `dx` from the Zygote format to the ChainRules differential types.
This is similar to `wrap_chainrules_input(dx)`, but because it gets `primal::T`,
it can turn `NamedTuple`s into `Tangent{T}(...)` not `Tangent{Any}(...)`.
"""
zygote2differential(x, primal) = z2d(x, primal)
zygote2differential(::Nothing, ::Any) = NoTangent()
zygote2differential(t::Tuple, primal::Tuple) = map(z2d, t, primal)
zygote2differential(t::Tuple, primal) = (@warn "primal should be a tuple, not $primal"; return t)

z2d(::Nothing, ::Any) = NoTangent()
z2d(::Tuple{Vararg{Nothing}}, ::Tuple) = NoTangent()  # collapse all-zero case
z2d(dx, ::Any) = dx
z2d(dx::AbstractArray{<:Number}, primal::AbstractArray) = dx
z2d(dx::AbstractArray{<:AbstractArray{<:Number}}, primal::AbstractArray) = dx
z2d(dx::AbstractArray, primal::AbstractArray) = map(z2d, dx, primal)
#=
# As an optimisation, we can convert by `reinterpret` for bitstypes, e.g. arrays of tuples of numbers
function z2d(dx::AbstractArray{S}, primal::AbstractArray{P}) where {S,P}
  if isbitstype(S)
    T = wrap_chainrules_input(S)
    reinterpret(Tangent{P,T}, dx)
  else
    map(z2d, dx, primal)
  end
end
=#

# Note: this should never be hit if we are converting things right, but it seems to be
# happening in the wild for sufficiently weird functions/types.
# This fixes most (all?) cases, but it would be good to find what we miss.
z2d(x::Union{AbstractZero, Tangent}, ::Any) = return x

function z2d(delta::Tuple, primal::Tuple)
  backing = map(z2d, delta, primal)
  if backing isa Tuple{Vararg{AbstractZero}}
    return NoTangent()  # collapse all-zero case
  else
    return canonicalize(Tangent{typeof(primal), typeof(backing)}(backing))
  end
end

# Dict handling in Zygote is a mess... should this become a  `Tangent{Dict,Dict}` ?
# Right now it uses a NamedTuple but not for fields of the AbstractDict struct
z2d(dx::NamedTuple, primal::AbstractDict) = dx

function z2d(delta::NamedTuple, primal::T) where T  # arbitrart struct
  fnames = fieldnames(T)
  deltas = map(n -> get(delta, n, nothing), fnames)
  primals = map(n -> getfield(primal, n), fnames)
  inner = map(z2d, deltas, primals)  # recurse into fields
    if inner isa Tuple{Vararg{AbstractZero}}
    return NoTangent()  # collapse all-zero case
  else
    backing = NamedTuple{fnames}(inner)
    return canonicalize(Tangent{T, typeof(backing)}(backing))
  end
end

# Dict case matches signature for ambiguity reasons:
z2d(dx::NamedTuple{L,S}, primal::AbstractDict) where {L,S<:Tuple{Vararg{Union{Number,Nothing}}}} = dx
# On Julia <= 1.6, this fixes easy cases which do not require recursion into fields, e.g.
# @inferred Zygote.z2d((re=1, im=nothing), 3.0+im)
@generated function z2d(delta::NamedTuple{L,S}, primal::T) where {L,S<:Tuple{Vararg{Union{Number,Nothing}}}, T}
  fnames = fieldnames(T)
  deltas = map(fnames) do n
    i = findfirst(isequal(n), L)
    if i == nothing || S.parameters[i] == Nothing
      :(NoTangent())
    else
      :(delta.$n)
    end
  end
  if all(d -> d == :(NoTangent()), deltas)
    return :(NoTangent())  # collapse all-zero case
  else
    return quote
      backing = NamedTuple{$fnames}(($(deltas...),))
      Tangent{$T, typeof(backing)}(backing)
    end
  end
end

z2d(dx::Ref, primal) = z2d(dx[], primal)  # mutable structs


text:
# we generate code in this module, so precompile where possible
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module Geodesy

using FixedSizeArrays
using Proj4

import Base.show, Base.call, Base.convert, Base.+, Base.-, Base.isnan, Base.getindex

# decalare the position supertype here
abstract  GeodesyType

for f in ["datums", "ellipsoids", "srid", "geodesy_types", "known_srids", "utm", "math_funcs", "type_methods", "transform", "external",
          "ITRF_transformations",
          "bounds", "vicenty", "distance",                                                          # do these belong here?
          "Proj4_Types/proj4_types", "Proj4_Types/projections", "Proj4_Types/proj4_transforms"]     # these definitely don't
    include("$f.jl")
end

export

    #
    # Position types supported by Geodesy directly
    #

    # World coordinate systems
    LLA,
    LL,
    ECEF,

    # Local coordinate systems
    ENU,

    #
    # Position types supported by Proj4
    #
    CRS,
    CCRS_Geoid,

    #
    # convenience for using points with a set datum
    #
    LLA_WGS84,   # typealias for LLA{WGS84}
    ECEF_WGS84,  # typealias for ECEF{WGS84}

    # no datum info
    LL_NULL,
    LLA_NULL,
    ECEF_NULL,
    CRS_NULL,
    CCRS_NULL,
    ENU_NULL,

    # Other types
    # Bounds,  # I dont want to export something named Bounds, maybe rename to Geobounds?

    #
    # export some datums
    #
    AbstractDatum,
    UnknownDatum,
    WGS84,
    GRS80,
    GDA94,
    NZGD2000,

    # datum functions
    get_datum,
    get_datums,    # get a list of all known datums

    #
    # transformation functions
    #
    geotransform,           # transform a point
    geotransform_vector,    # transform a vector of points
    geotransform_params,    # parameters for a transformation

    #
    # SRID related
    #
    AbstractSRID,
    UnknownSRID,
    SRID,
    get_srid,

    #
    # Geoid related
    #
    AbstractGeoid,
    UnknownGeoid,
    get_geoid,

    # utm related
    utm_zone,

    #
    # accessors
    #
    getX,       # should I really export this?
    getY,
    getZ,
    get_lat,
    get_lon,
    get_alt,
    get_east,
    get_north,
    get_up,

    # Methods
    # center,   # dont want to export names this generic
    # distance, # dont want to export names this generic

    # convert to / from geodesy types
    geodify,
    ungeodify,

    decimal2dms,
    dms2decimal

    #=
    inBounds,
    haversine_distance,
    boundaryPoint,
    onBounds
    =#

end # module Geodesy


text:
@doc raw"""
    HyperbolicDiffusionEquations1D

The linear hyperbolic diffusion equations in one space dimension.
A description of this system can be found in Sec. 2.5 of the book
- Masatsuka (2013)
  I Do Like CFD, Too: Vol 1.
  Freely available at [http://www.cfdbooks.com/](http://www.cfdbooks.com/)
Further analysis can be found in the paper
- Nishikawa (2007)
  A first-order system approach for diffusion equation. I: Second-order residual-distribution
  schemes
  [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029)
"""
struct HyperbolicDiffusionEquations1D{RealT<:Real} <: AbstractHyperbolicDiffusionEquations{1, 2}
  Lr::RealT     # reference length scale
  inv_Tr::RealT # inverse of the reference time scale
  nu::RealT     # diffusion constant
end

function HyperbolicDiffusionEquations1D(; nu=1.0, Lr=inv(2pi))
  Tr = Lr^2 / nu
  HyperbolicDiffusionEquations1D(promote(Lr, inv(Tr), nu)...)
end

get_name(::HyperbolicDiffusionEquations1D) = "HyperbolicDiffusionEquations1D"
varnames(::typeof(cons2cons), ::HyperbolicDiffusionEquations1D) = ("phi", "q1")
varnames(::typeof(cons2prim), ::HyperbolicDiffusionEquations1D) = ("phi", "q1")
default_analysis_errors(::HyperbolicDiffusionEquations1D) = (:l2_error, :linf_error, :residual)

@inline function residual_steady_state(du, ::HyperbolicDiffusionEquations1D)
  abs(du[1])
end

"""
    initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)

A non-priodic smooth initial condition. Can be used for convergence tests in combination with
[`source_terms_poisson_nonperiodic`](@ref) and [`boundary_condition_poisson_nonperiodic`](@ref).
!!! note
    The solution is periodic but the initial guess is not.
"""
function initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)
  # elliptic equation: -νΔϕ = f
  # Taken from Section 6.1 of Nishikawa https://doi.org/10.1016/j.jcp.2007.07.029
  if t == 0.0
    # initial "guess" of the solution and its derivative
    phi = x[1]^2 - x[1]
    q1  = 2*x[1] - 1
  else
    phi = sinpi(x[1])      # ϕ
    q1  = pi * cospi(x[1]) # ϕ_x
  end
  return SVector(phi, q1)
end

"""
    source_terms_poisson_nonperiodic(u, x, t,
                                     equations::HyperbolicDiffusionEquations1D)

Source terms that include the forcing function `f(x)` and right hand side for the hyperbolic
diffusion system that is used with [`initial_condition_poisson_nonperiodic`](@ref) and
[`boundary_condition_poisson_nonperiodic`](@ref).
"""
@inline function source_terms_poisson_nonperiodic(u, x, t,
                                                  equations::HyperbolicDiffusionEquations1D)
  # elliptic equation: -νΔϕ = f
  # analytical solution: ϕ = sin(πx) and f = π^2sin(πx)
  @unpack inv_Tr = equations

  dphi = pi^2 * sinpi(x[1])
  dq1  = -inv_Tr * u[2]

  return SVector(dphi, dq1)
end

"""
    boundary_condition_poisson_nonperiodic(u_inner, orientation, direction, x, t,
                                        surface_flux_function,
                                        equations::HyperbolicDiffusionEquations1D)

Boundary conditions used for convergence tests in combination with
[`initial_condition_poisson_nonperiodic`](@ref) and [`source_terms_poisson_nonperiodic`](@ref).
"""
function boundary_condition_poisson_nonperiodic(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equations::HyperbolicDiffusionEquations1D)
  # elliptic equation: -νΔϕ = f
  phi = sinpi(x[1])      # ϕ
  q1  = pi * cospi(x[1]) # ϕ_x
  u_boundary = SVector(phi, q1)

  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)

A non-priodic harmonic function used in combination with
[`source_terms_poisson_nonperiodic`](@ref) and [`boundary_condition_poisson_nonperiodic`](@ref).
!!! note
    The only harmonic functions in 1D have the form phi(x) = A + Bx
"""
function initial_condition_harmonic_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D)
  # elliptic equation: -νΔϕ = f
  if t == 0.0
    phi = 5.0
    q1  = 0.0
  else
    A = 3
    B = exp(1)
    phi = A + B * x[1]
    q1  = B
  end
  return SVector(phi, q1)
end

"""
    source_terms_harmonic(u, x, t, equations::HyperbolicDiffusionEquations1D)

Source term that only includes the forcing from the hyperbolic diffusion system
used with [`initial_condition_harmonic_nonperiodic`](@ref) and
[`boundary_condition_harmonic_nonperiodic`](@ref).
"""
@inline function source_terms_harmonic(u, x, t, equations::HyperbolicDiffusionEquations1D)
  # harmonic solution of the form ϕ = A + B * x, so f = 0
  @unpack inv_Tr = equations

  dq1 = -inv_Tr * u[2]

  return SVector(zero(dq1), dq1)
end

"""
    boundary_condition_harmonic_nonperiodic(u_inner, orientation, direction, x, t,
                                        surface_flux_function,
                                        equations::HyperbolicDiffusionEquations1D)

Boundary conditions used for the harmonic function test case
[`initial_condition_harmonic_nonperiodic`](@ref) and [`source_terms_harmonic`](@ref).
"""
function boundary_condition_harmonic_nonperiodic(u_inner, orientation, direction, x, t,
                                                 surface_flux_function,
                                                 equations::HyperbolicDiffusionEquations1D)
  u_boundary =  initial_condition_harmonic_nonperiodic(x, t, equations)

  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations1D)

Setup used for convergence tests of the Euler equations with self-gravity used in
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
in combination with [`source_terms_harmonic`](@ref).
"""
function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations1D)

  # Determine phi_x
  G = 1.0           # gravitational constant
  C = -4.0 * G / pi # -4 * G / ndims * pi
  A = 0.1           # perturbation coefficient must match Euler setup
  rho1 = A * sinpi(x[1] - t)
  # intialize with ansatz of gravity potential
  phi = C * rho1
  q1  = C * A * pi * cospi(x[1] - t) # = gravity acceleration in x-direction

  return SVector(phi, q1)
end

# Calculate 1D flux in for a single point
@inline function flux(u, orientation, equations::HyperbolicDiffusionEquations1D)
  phi, q1 = u
  @unpack inv_Tr = equations

  # Ignore orientation since it is always "1" in 1D
  f1 = -equations.nu * q1
  f2 = -phi * inv_Tr

  return SVector(f1, f2)
end

@inline function flux_lax_friedrichs(u_ll, u_rr, orientation, equations::HyperbolicDiffusionEquations1D)
  # Obtain left and right fluxes
  f_ll = flux(u_ll, orientation, equations)
  f_rr = flux(u_rr, orientation, equations)

  λ_max = sqrt(equations.nu * equations.inv_Tr)

  return 0.5 * (f_ll + f_rr - λ_max * (u_rr - u_ll))
end

@inline have_constant_speed(::HyperbolicDiffusionEquations1D) = Val(true)

@inline function max_abs_speeds(eq::HyperbolicDiffusionEquations1D)
  return sqrt(eq.nu * eq.inv_Tr)
end

# Convert conservative variables to primitive
@inline cons2prim(u, equations::HyperbolicDiffusionEquations1D) = u

# Convert conservative variables to entropy found in I Do Like CFD, Too, Vol. 1
@inline function cons2entropy(u, equations::HyperbolicDiffusionEquations1D)
  phi, q1 = u

  w1 = phi
  w2 = equations.Lr^2 * q1

  return SVector(w1, w2)
end

# Calculate entropy for a conservative state `u` (here: same as total energy)
@inline entropy(u, equations::HyperbolicDiffusionEquations1D) = energy_total(u, equations)

# Calculate total energy for a conservative state `u`
@inline function energy_total(u, equations::HyperbolicDiffusionEquations1D)
  # energy function as found in equations (2.5.12) in the book "I Do Like CFD, Vol. 1"
  phi, q1 = u
  return 0.5 * (phi^2 + equations.Lr^2 * q1^2)
end


text:
using Yao, LinearAlgebra, YaoBlocks
import Random.shuffle

Ugate(nbit::Int, i::Int) = put(nbit, i=>chain(Rz(0), Rx(0), Rz(0)));
Ugate(i::Int) = put(i=>chain(Rz(0), Rx(0), Rz(0)));

struct InitStateMatrix
    m1::Matrix{ComplexF64}
    m2::Matrix{ComplexF64}
end

function build_state(mats::InitStateMatrix)
    # return zero_state(2) |> chain(2, put(1=>chain(Rz(1), Rx(2))), put(2=>chain(Rz(1), Rx(2))))
    return zero_state(3) |> chain(3, put(3, (2)=>matblock(mats.m1)), put(3, (3)=>matblock(mats.m2)))
end

# @show probs(build_state())

function overlap(mats::InitStateMatrix)
    return dot(zero_state(2) |> chain(2, put(2, (1)=>matblock(mats.m1)), put(2, (2)=>matblock(mats.m2)), control(1, 2=>X), put(1=>H)) |> probs,[1, 1, 1, -1])
end

function cost(trainX, trainY, circuit, c)
    actual_out = [dot(i |> build_state |> circuit |> probs, c) for i in trainX]
    return sum((trainY-actual_out).^2)/length(trainY)
end

function generate_gates_set(n::Int64)
    rtn = []
    for i in 1:n
        for j in 1:n
            if i != j
                push!(rtn, control(n, i, j=>X))
            end
        end
    end

    for i in 1:n
        push!(rtn, Ugate(n, i))
    end

    return rtn
end

# 检测生成的circuit是否合法
function validate_circuit(circuit)
    old = circuit[1]
    for i in 2:length(circuit)
        now = circuit[i]
        if isa(old, ControlBlock) && old == now
            return false
        end
        if (!isa(old, ControlBlock)) && (!isa(now, ControlBlock)) && occupied_locs(old) == occupied_locs(now)
            return false
        end
        old = now
    end

    if length(circuit) > 3 && length(param_indexes(circuit)) == 0
        return false
    end

    return true
end

function param_indexes(circuit)
    rtn = []
    for i in 1:length(circuit)
        if nparameters(circuit[i]) > 0
            push!(rtn, i)
        end
    end
    return rtn
end

function generate_init_circuit(n::Int64, d::Int64)
    gate_set = generate_gates_set(n)
    rtn = chain(n)
    while length(rtn) < d
        push!(rtn, rand(gate_set))
        if !validate_circuit(rtn)
            pop!(rtn)
        end
    end
    params = rand(nparameters(rtn)) * 2π .- π
    dispatch!(rtn, params)

    return rtn
end

function random_modify(circuit, circuit_indexes, gate_set)
    circuit2 = deepcopy(circuit)
    index = rand(circuit_indexes)
    # @show index
    circuit2[index] = rand(gate_set)
    if validate_circuit(circuit2)
        params = rand(nparameters(circuit2[index])) * 2π .- π
        dispatch!(circuit2[index], params)
        return circuit2
    end
    return random_modify(circuit, circuit_indexes, gate_set)
end

function gradient(circuit, circuit_indexes, trainX, trainY, delta)
    n = nqubits(circuit)
    grad = zeros(Float64, nparameters(circuit))

    count = 1
    for i in 1:length(circuit)
        params = parameters(circuit[i])
        for k in 1:nparameters(circuit[i])
            if i in circuit_indexes
                params[k] += delta
                dispatch!(circuit[i], params)
                cost_pos = cost(trainX, trainY, circuit, c)

                params[k] -= 2 * delta
                dispatch!(circuit[i], params)
                cost_neg = cost(trainX, trainY, circuit, c)

                params[k] += delta
                dispatch!(circuit[i], params)

                # @show (cost_pos - cost_neg)

                grad[count] = (cost_pos - cost_neg)/(2 * delta)
            end
            count += 1
        end
    end
    return grad
end

using Flux.Optimise

c = [1, -1, 1, -1, 1, -1, 1, -1];
const cost_threshold = 2e-7
# circuit = build_circuit()
trainX = [InitStateMatrix(rand_unitary(2), rand_unitary(2)) for _ in 1:10]
testX = [InitStateMatrix(rand_unitary(2), rand_unitary(2)) for _ in 1:10]
trainY = map(overlap, trainX)
testY = map(overlap, testX)
opt = Descent(.9)
# opt = ADAM(0.03, (0.9, 0.999))

function train(circuit, circuit_indexes, trainX, trainY, testX, testY, cost_threshold)
    history = Float64[]
    history2 = Float64[]
    last_cost = now_cost = 0
    for i in 1:50
        last_cost = now_cost
        now_cost = cost(testX, testY, circuit, c)
        push!(history, cost(trainX, trainY, circuit, c))
        push!(history2, now_cost)
        ps = parameters(circuit)
        # @show (i, ps, now_cost)
        if (abs(now_cost - last_cost) < now_cost * .005) || now_cost < cost_threshold
            # @show now_cost
            break
        end

        # prevent local-minimal
        for index in shuffle(circuit_indexes)
            Optimise.update!(opt, ps, gradient(circuit, [ index ], trainX, trainY, .001))
        end
        popdispatch!(circuit, ps)
        # if i == 200
        #     println("Warning: failed to converge")
        # end
    end
    return (now_cost, history, history2)
end

test_circuit = generate_init_circuit(3, 9)
(tmp, history, history2) = train(test_circuit, param_indexes(test_circuit), trainX, trainY, testX, testY, cost_threshold)
# (tmp, history, history2) = train(test, [ rand(param_indexes(test)) ], trainX, trainY, testX, testY, cost_threshold)
using Plots
function plot_history(history)
    pyplot() # Choose a backend
    fig1 = plot([history, history2]; legend=nothing)
    title!("training history")
    xlabel!("steps"); ylabel!("cost")
end
plot_history(history)


text:
##########
# Others #
##########

# VarInfo to Sample
@inline Sample(vi::VarInfo) = begin
  value = Dict{Symbol, Any}() # value is named here because of Sample has a field called value
  for vn in keys(vi)
    value[sym(vn)] = realpart(vi[vn])
  end

  # NOTE: do we need to check if lp is 0?
  value[:lp] = realpart(getlogp(vi))

  if ~isempty(vi.pred)
    for sym in keys(vi.pred)
      # if ~haskey(sample.value, sym)
        value[sym] = vi.pred[sym]
      # end
    end
    # TODO: check why 1. 2. cause errors
    # TODO: which one is faster?
    # 1. Using empty!
    # empty!(vi.pred)
    # 2. Reassign an enmtpy dict
    # vi.pred = Dict{Symbol,Any}()
    # 3. Do nothing?
  end

  Sample(0.0, value)
end

# VarInfo, combined with spl.info, to Sample
@inline Sample(vi::VarInfo, spl::Sampler) = begin
  s = Sample(vi)

  if haskey(spl.info, :ϵ)
    s.value[:epsilon] = spl.info[:ϵ][end]
  end

  if haskey(spl.info, :lf_num)
    s.value[:lf_num] = spl.info[:lf_num][end]
  end

  s
end


text:
"""
    extractnotebook(notebook)
Reads a Pluto notebook file, extracts the code
and returns a dictionary with the code in organized form.
The output is a dictionary containing
* `codes`        - The cell code of each running cell;
* `notebookname` - Name of the notebook;
* `notebookdir`  - Directory location of the notebook file;
* `contents`     - Only the Julia code in the cell;
* `outputtag`    - Whether the output is hidden or showing (read the `tagcelloutput()` function);
* `celltype`     - Whether cell contains code or markdown;
* `order`        - Order that the cells are displayed in the notebook;
* `view`         - Whether the code is hidden or showing (the "eye" icon in the notebook);

e.g. `extractnotebook("./mynotebook.jl")`
notebookdata = Dict(:codes => codes, :notebookname => notebookname, :notebookdir => notebookdir,
                        :contents => contents, :outputtag=>outputtag,
                        :celltype => celltype,:order=> order,:view=>view)
"""
function extractnotebook(notebook)
    s = read(notebook, String)
    cells = split(s, "# ╔═╡ ");
    # The first cell and the final 3 are not used
    codes, contents, outputtag, celltype = [],[],[],[]
    for cell in cells[2:end-3]
        push!(codes, cell[1:36])
        push!(contents, cell[38:end])
        push!(outputtag, endswith(rstrip(cell),";") ? "hideoutput" : "showoutput")
        #= if length(cell) < 42 =#
        #=     push!(celltype, "code") =#
        #= else =#
        push!(celltype, startswith(strip(cell[38:end]),"md\"") ? "markdown" : "code")
    end

    # Get order and view type
    r = either("# ╠","# ╟")
    sortedcells = split(cells[end],Regex(r))
    sortedcodes = [cell[4:39] for cell in sortedcells[2:end-2]]
    order = [findfirst(isequal(scode),codes) for scode in  sortedcodes[1:end]]
    view  = [occursin("═",c) ? "showcode" : "hidecode" for c in sortedcells[2:end-2]]
    # Matching running order
    view  = view[[findfirst(isequal(scode),sortedcodes) for scode in  codes]]

    # inferring the notebook name
    # base on the notebook file path.
    notebookname = split(notebook,"/")[end]
    if endswith(notebookname,".jl")
        # removes the ".jl" in the end
        notebookname = notebookname[1:end-3]
    end
    notebookdir = dirname(notebook) == "" ? "./" : dirname(notebook)

    notebookdata = Dict(:codes => codes, :notebookname => notebookname, :notebookdir => notebookdir,
                        :contents => contents, :outputtag=>outputtag,
                        :celltype => celltype,:order=> order,:view=>view)
    return notebookdata
end

"""
    collectoutputs(notebookdata, path)
Runs the Pluto notebook and saves the outputs of
each cell into a list.
"""
function collectoutputs(notebookdata, path)
    figureindex = Dict(:i => 0)
    runpath = pwd()
    cd(notebookdata[:notebookdir])
    outputs = []
    for (i, content) ∈ enumerate(notebookdata[:contents])
        if notebookdata[:celltype][i] == "code"
            if startswith(lstrip(content),"begin") && endswith(rstrip(content),"end")
                ex = :($(Meta.parse(strip(content))))
            else
                ex = :($(Meta.parse("begin\n"*content*"\n end")))
            end

            s = string(ex.args[end])
            if contains(s, Regex(either("PlutoUI.LocalResource","LocalResource")))
                if findfirst(Regex(look_for(one_or_more(ANY),after="(\"",before="\")")),s) === nothing
                    Core.eval(Runner, ex)
                    pathvariable = s[findfirst(Regex(look_for(one_or_more(ANY),after="(",before=")")),s)]
                    imagepath = Core.eval(Runner,Meta.parse(pathvariable))
                else
                    imagepath = s[findfirst(Regex(look_for(one_or_more(ANY),after="(\"",before="\")")),s)]
                end
                push!(outputs,(:image,pwd()*"/"*imagepath))
            else
                io = IOBuffer();
                cd(runpath)

                Base.invokelatest(show,
                    IOContext(io, :limit => true),"text/plain",
                    dispatch_output(Core.eval(Runner,ex), notebookdata[:notebookname], path, figureindex));
                cd(notebookdata[:notebookdir])

                celloutput = String(take!(io))
                if celloutput == "nothing"
                    push!(outputs,(:nothing, ""))
                elseif startswith(celloutput, "Plot{Plots.")
                    push!(outputs,
                    (:plot,notebookdata[:notebookname]*"_"*"figure"*string(figureindex[:i])*".png"))
                elseif startswith(celloutput, "FigureAxisPlot()")
                    push!(outputs,
                    (:plot,notebookdata[:notebookname]*"_"*"figure"*string(figureindex[:i])*".pdf"))
                else
                    push!(outputs,(:text, celloutput))
                end
            end
        else
            push!(outputs,nothing)
        end
    end
    cd(runpath)
    return outputs
end

function dispatch_output(command_eval::Plots.Plot, notebookname, path, figureindex)
    figureindex[:i]+=1
    Plots.savefig(command_eval,path*"/figures/"*notebookname*"_"*"figure"*string(figureindex[:i])*".png")
    return command_eval
end

function dispatch_output(command_eval::Makie.FigureAxisPlot, notebookname, path, figureindex)
    figureindex[:i]+=1
    CairoMakie.save(path*"/figures/"*notebookname*"_"*"figure"*string(figureindex[:i])*".pdf", command_eval)
    return command_eval
end

"""
    plutotolatex(notebookname, targetdir="./build_latex"; template=:book, fontpath=nothing)
Function to convert Pluto notebooks. The arguments are the same as the ones in
`notebooktolatex`.
"""
function plutotolatex(notebookname, targetdir="./build_latex"; template=:book, fontpath=nothing)

    createproject(targetdir, template, fontpath)
    nb = extractnotebook(notebookname)
    texfile = read(targetdir*"/main.tex", String)
    lineinsert = 1
    for (i,line) in enumerate(split(texfile, "\n"))
            lineinsert = i
            if startswith(line, "% INCLUDE NOTEBOOKS")
            break
        end
    end

    if !occursin("\\include{./notebooks/"*nb[:notebookname]*"}",read(targetdir*"/main.tex", String))

        insertlinebelow(targetdir*"/main.tex",
            "\\include{./notebooks/"*nb[:notebookname]*"}", lineinsert)
    end

    outputs = collectoutputs(nb,targetdir);
    notebook = targetdir*"/notebooks/"*nb[:notebookname]*".tex"
    open(notebook, "w") do f
        write(f,"\\newpage\n")
        for i in nb[:order]
            if nb[:celltype][i] == "markdown"
                if startswith(strip(nb[:contents][i]), "md\"\"\"")
                    parsed = markdowntolatex(strip(nb[:contents][i])[7:end-3],
                        targetdir, nb[:notebookdir], template=template)*"\n\n"
                elseif startswith(strip(nb[:contents][i]), "md\"")
                    parsed = markdowntolatex(strip(nb[:contents][i])[4:end-1],
                        targetdir, nb[:notebookdir], template=template)*"\n\n"
                else
                    throw(DomainError("Markdown cell must start with either md\"\"\" or md\"."))
                end
                write(f,parsed)
            elseif nb[:celltype][i] == "code" && nb[:view][i] == "showcode"
                write(f,"\n\\begin{lstlisting}[language=JuliaLocal, style=julia]\n")
                write(f, strip(nb[:contents][i]))
                write(f,"\n\\end{lstlisting}\n")
            end
            if nb[:celltype][i] == "code" && nb[:outputtag][i] == "showoutput"
                if outputs[i][1] == :text
                    write(f,"\n\\begin{verbatim}\n")
                    write(f, outputs[i][2])
                    write(f,"\n\\end{verbatim}\n")
                elseif outputs[i][1] == :plot
                    write(f,"\n\\begin{figure}[H]\n")
                    write(f,"\t\\centering\n")
                    write(f,"\t\\includegraphics[width=0.8\\textwidth]{./figures/"*outputs[i][2]*"}\n")
                    write(f,"\t\\label{fig:"*outputs[i][2]*"}\n")
                    write(f,"\n\\end{figure}\n")
                elseif outputs[i][1] == :image
                    figurename = ""
                    if outputs[i][2][end-3:end] == ".svg"
                        figuresvg = outputs[i][2]
                        figurename *= basename(figuresvg[1:end-4])
                        figurepdf = targetdir*"/figures/"*figurename*".pdf"
                        rsvg_convert() do cmd
                           run(`$cmd $figuresvg -f pdf -o $figurepdf`)
                        end
                    else
                        figurename *= basename(outputs[i][2])
                        cp(outputs[i][2],targetdir*"/figures/"*figurename, force=true)
                    end
                    write(f,"\n\\begin{figure}[H]\n")
                    write(f,"\t\\centering\n")
                    write(f,"\t\\includegraphics[width=0.8\\textwidth]{./figures/"*figurename*"}\n")
                    #= write(f,"\t\\includegraphics[width=0.8\\textwidth]{"*outputs[i][2]*"}\n") =#
                    write(f,"\t\\label{fig:"*outputs[i][2]*"}\n")
                    write(f,"\n\\end{figure}\n")
                end
            end
        end
    end
end


text:
function chooseproc()
    nextidx = 0
    while true
        p = workers()[(nextidx) + 1]
        nextidx += 1
        if nextidx >= nworkers()
            nextidx = 0
        end
        produce(p)
    end
end


text:
#=
Script to generate QASM for Quantum Fourier Transform circuits
=#

using PyCall
using ArgParse
using PicoQuant

"""
    parse_commandline()

Parse command line options and return argument dictionary
"""
function parse_commandline(ARGS)
    s = ArgParseSettings()
        @add_arg_table! s begin
            "--number-qubits", "-n"
                help = "Size of circuit"
                arg_type = Int
                required = true
            "--output", "-o"
                help = "Output file to write tensor network graph to"
                arg_type = String
                default = ""
            "--latex-circuit"
                help = "Generate latex for circuit diagram"
                action = :store_true
            "--text-circuit"
                help = "Generate ASCII for circuit diagram"
                action = :store_true
        end
        return parse_args(ARGS, s)
end

function main(ARGS)
    parsed_args = parse_commandline(ARGS)

    qubits = parsed_args["number-qubits"]

    circuit = create_qft_circuit(qubits)

    if parsed_args["output"] == ""
        filename = "qft_$(qubits).qasm"
    else
        filename = parsed_args["output"]
    end

    open(filename, "w") do io
        write(io, circuit.qasm())
    end

    if parsed_args["latex-circuit"]
        latex_filename = "$(splitext(filename)[1]).tex"
        circuit.draw(output="latex_source", filename=latex_filename)
    end

    if parsed_args["text-circuit"]
        text_filename = "$(splitext(filename)[1]).txt"
        circuit.draw(output="text", filename=text_filename)
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text:
# Fit mixture model using MOSEK interior-point solver, called by the
# KWDual function from the REBayes package.
function REBayes(L)
  @rput L;
  R"library(REBayes);
    n      <- nrow(L)
    k      <- ncol(L)
    timing <- system.time(out <- KWDual(L,rep(1,k),rep(1,n)/n));
    timing <- timing[3];
    x      <- out$f;
    x      <- x/sum(x);"
  @rget x;
  @rget timing;
  return x, timing
end

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper