DominantCycle
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © lastguru

//@version=5
// @description Collection of Dominant Cycle estimators. Length adaptation used in the Adaptive Moving Averages and the Adaptive Oscillators try to follow price movements and accelerate/decelerate accordingly (usually quite rapidly with a huge range). Cycle estimators, on the other hand, try to measure the cycle period of the current market, which does not reflect price movement or the rate of change (the rate of change may also differ depending on the cycle phase, but the cycle period itself usually changes slowly). This collection may become encyclopaedic, so if you have any working cycle estimator, drop me a line in the comments below. Suggestions are welcome. Currently included estimators are based on the work of John F. Ehlers
library("DominantCycle", true)
import lastguru/CommonFilters/1 as cf
import lastguru/NormalizedOscillators/1 as co
import lastguru/LengthAdaptation/3 as la

// ----------
// Dominant Cycle Estimators
// ----------

// @function MESA Adaptation - MAMA Cycle
// @param src Series to use
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @returns Calculated period
// Based on MESA Adaptive Moving Average by John F. Ehlers
// Performs Hilbert Transform Homodyne Discriminator cycle measurement
// Unlike MAMA Alpha function (in LengthAdaptation library), this does not compute phase rate of change
// Introduced in the September 2001 issue of Stocks and Commodities
// Inspired by the @everget implementation: https://www.tradingview.com/script/aaWzn9bK-Ehlers-MESA-Adaptive-Moving-Averages-MAMA-FAMA/
// Inspired by the @anoojpatel implementation: https://www.tradingview.com/script/uUHctAeu-Missile-RSI-RSI-of-momentum-w-Dominant-Cycle-length-Fisher/
export mamaPeriod(float src, int dynLow, int dynHigh) =>
    period = 0.0
    C1 = 0.0962
    C2 = 0.5769
    C3 = 0.075 * nz(period[1]) + 0.54
    
    smooth = (4 * src + 3 * nz(src[1]) + 2 * nz(src[2]) + nz(src[3])) / 10
    detrend = C3 * (C1 * smooth + C2 * nz(smooth[2]) - C2 * nz(smooth[4]) - C1 * nz(smooth[6]))

    // Compute InPhase and Quadrature components
    Q1 = C3 * (C1 * detrend + C2 * nz(detrend[2]) - C2 * nz(detrend[4]) - C1 * nz(detrend[6]))
    I1 = nz(detrend[3])

    // Advance Phase of I1 and Q1 by 90 degrees
    jI = C3 * (C1 * I1 + C2 * nz(I1[2]) - C2 * nz(I1[4]) - C1 * nz(I1[6]))
    jQ = C3 * (C1 * Q1 + C2 * nz(Q1[2]) - C2 * nz(Q1[4]) - C1 * nz(Q1[6]))

    // Phasor addition for 3 bar averaging
    I2 = I1 - jQ
    Q2 = Q1 + jI

    // Smooth I and Q components before applying discriminator
    I2 := 0.2*I2 + 0.8*nz(I2[1])
    Q2 := 0.2*Q2 + 0.8*nz(Q2[1])

    // Extract Homodyne Discriminator
    Re = I2*nz(I2[1]) + Q2*nz(Q2[1])
    Im = I2*nz(Q2[1]) - Q2*nz(I2[1])
    Re  := 0.2*Re + 0.8*nz(Re[1])
    Im  := 0.2*Im + 0.8*nz(Im[1]) 

    period := (Re != 0 and Im != 0) ? 2 * math.pi / math.atan(Im/Re) : 0
    period := math.min(period, 1.5 * nz(period[1], period))
    period := math.max(period, (2/3) * nz(period[1], period))
	period := math.min(math.max(period, dynLow), dynHigh)
    period := period*0.2 + nz(period[1])*0.8
    period

// @function Pearson Autocorrelation
// @param src Series to use
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @param preHP Use High Pass prefilter (default)
// @param preSS Use Super Smoother prefilter (default)
// @param preHP Use Hann Windowing prefilter
// @returns Calculated period
// Based on Pearson Autocorrelation Periodogram by John F. Ehlers
// Introduced in the September 2016 issue of Stocks and Commodities
// Inspired by the @blackcat1402 implementation: https://www.tradingview.com/script/9CbjlKjz-blackcat-L2-Ehlers-Adaptive-CCI-2013/
// Inspired by the @rumpypumpydumpy implementation: https://www.tradingview.com/script/RV0MQvIB-Ehler-s-Autocorrelation-Periodogram-RSI-MFI/
// Corrected many errors, and made small speed optimizations, so this could be the best implementation to date (still slow, though, so may revisit in future)
// High Pass and Super Smoother prefilters are used in the original implementation
export paPeriod (float src, int dynLow, int dynHigh, bool preHP = true, bool preSS = true, bool preHW = false) =>
    avglen = 3
    var corr = array.new_float(dynHigh+1, initial_value=0)
    var R = array.new_float(dynHigh+1, initial_value=0)
    var pwr = array.new_float(dynHigh+1, initial_value=0)

    hp = preHP ? cf.hp(src, dynHigh) : src
    pfilt = preSS ? cf.supers2(hp, dynLow) : hp
    filt = preHW ? cf.hannFIR(pfilt, dynLow) : pfilt

    //Pearson correlation for each value of lag
    for lag = 0 to dynHigh
    	//Set the averaging length as M 
    	m = avglen == 0 ? lag : avglen
    	Sx = 0.0
    	Sy = 0.0
    	Sxx = 0.0
    	Sxy = 0.0
    	Syy = 0.0
        
    	//Advance samples of both data streams and sum Pearson components
    	for count = 0 to m-1
    	    x = filt[count]
    	    y = nz(filt[lag + count])
    	    Sx := Sx + x
    	    Sy := Sy + y
    	    Sxx := Sxx + x * x
    	    Sxy := Sxy + x * y
    	    Syy := Syy + y * y
        
    	//Compute correlation for each value of lag
    	array.set(corr, lag, (m*Sxx - Sx*Sx)*(m*Syy - Sy*Sy) > 0 ? (m*Sxy - Sx*Sy) / math.sqrt((m*Sxx - Sx*Sx) * (m*Syy - Sy*Sy)) : 0)
    
    //Compute the Fourier Transform for each Correlation
    for period = dynLow to dynHigh
        cosinePart = 0.0
        sinePart = 0.0
    	for n = avglen to dynHigh
    		cosinePart := cosinePart + array.get(corr, n) * math.cos(2 * math.pi * n / period)
    		sinePart := sinePart + array.get(corr, n) * math.sin(2 * math.pi * n / period)
    	sqSum = math.pow(cosinePart, 2) + math.pow(sinePart, 2)
    	array.set(R, period, 0.2 * math.pow(sqSum, 2) + 0.8 * nz(array.get(R, period), math.pow(sqSum, 2)))

    //Find Maximum Power Level for Normalization
    maxPwr = 0.0
    for period = dynLow to dynHigh
        maxPwr := math.max(array.get(R, period), maxPwr)
    
    for period = dynLow to dynHigh
    	array.set(pwr, period, array.get(R, period) / maxPwr)
    
    //Compute the dominant cycle using the CG of the spectrum
    peakPwr = 0.0
    for period = dynLow to dynHigh
        peakPwr := math.max(array.get(pwr, period), peakPwr)
    
    Spx = 0.0
    Sp = 0.0
    for period = dynLow to dynHigh
        if peakPwr >= 0.25 and array.get(pwr, period) >= 0.25
            Spx := Spx + period * array.get(pwr, period)
            Sp := Sp + array.get(pwr, period)
    
    domCycle = 0.0
    if Sp != 0
        domCycle := Spx / Sp
    if Sp < 0.25
        domCycle := nz(domCycle[1], domCycle)
    
    domCycle

// @function Discrete Fourier Transform
// @param src Series to use
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @param preHP Use High Pass prefilter (default)
// @param preSS Use Super Smoother prefilter (default)
// @param preHP Use Hann Windowing prefilter
// @returns Calculated period
// Based on Spectrum from Discrete Fourier Transform by John F. Ehlers
// Inspired by the @blackcat1402 implementation: https://www.tradingview.com/script/V6wCdSwQ-blackcat-L2-Ehlers-DFT-Adapted-RSI/
// High Pass, Super Smoother and Hann Windowing prefilters are used in the original implementation
export dftPeriod (float src, int dynLow, int dynHigh, bool preHP = true, bool preSS = true, bool preHW = false) =>
    var pwr = array.new_float(dynHigh+1, initial_value=0)

    hp = preHP ? cf.hp(src, dynHigh) : src
    pfilt = preSS ? cf.supers2(hp, dynLow) : hp
    filt = preHW ? cf.hannFIR(pfilt, dynLow) : pfilt

    //Compute the Fourier Transform
    for period = dynLow to dynHigh
        cosinePart = 0.0
        sinePart = 0.0
    	for n = 0 to dynHigh - 1
    		cosinePart := cosinePart + nz(filt[n]) * math.cos(2 * math.pi * n / period)
    		sinePart := sinePart + nz(filt[n]) * math.sin(2 * math.pi * n / period)
    	array.set(pwr, period, math.pow(cosinePart, 2) + math.pow(sinePart, 2))

    //Find Maximum Power Level for Normalization
    maxPwr = array.get(pwr, dynLow)
    for period = dynLow to dynHigh
        maxPwr := math.max(array.get(pwr, period), maxPwr)
    
    //Normalize Power Levels and Convert to Decibels
    for period = dynLow to dynHigh
    	array.set(pwr, period, array.get(pwr, period) / maxPwr)
    
    //Compute the dominant cycle using the CG of the spectrum
    Spx = 0.0
    Sp = 0.0
    for period = dynLow to dynHigh
        if array.get(pwr, period) >= 0.5
            Spx := Spx + period * array.get(pwr, period)
            Sp := Sp + array.get(pwr, period)
    
    domCycle = 0.0
    if Sp != 0
        domCycle := Spx / Sp

    domCycle

// @function Phase Accumulation
// @param src Series to use
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @param preHP Use High Pass prefilter (default)
// @param preSS Use Super Smoother prefilter (default)
// @param preHP Use Hamm Windowing prefilter
// @returns Calculated period
// Based on Dominant Cycle from Phase Accumulation by John F. Ehlers
// High Pass and Super Smoother prefilters are used in the original implementation
export phasePeriod (float src, int dynLow, int dynHigh, bool preHP = true, bool preSS = true, bool preHW = false) =>
    hp = preHP ? cf.hp(src, dynHigh) : src
    pfilt = preSS ? cf.supers2(hp, dynLow) : hp
    filt = preHW ? cf.hannFIR(pfilt, dynLow) : pfilt

    //Correlate with Cosine having a fixed period
	Sx = 0.0
	Sy = 0.0
	Sxx = 0.0
	Sxy = 0.0
	Syy = 0.0
    
	for count = 0 to dynHigh-1
	    x = nz(filt[count])
	    y = math.cos(2 * math.pi * count / dynHigh)
	    Sx := Sx + x
	    Sy := Sy + y
	    Sxx := Sxx + x * x
	    Sxy := Sxy + x * y
	    Syy := Syy + y * y
    
	real = (dynHigh * Sxx - Sx*Sx > 0) and (dynHigh * Syy - Sy*Sy > 0) ? (dynHigh * Sxy - Sx*Sy) / math.sqrt((dynHigh * Sxx - Sx*Sx) * (dynHigh * Syy - Sy*Sy)) : 0
	
    //Correlate with negative Sine having a fixed period
	Sx := 0.0
	Sy := 0.0
	Sxx := 0.0
	Sxy := 0.0
	Syy := 0.0
    
	for count = 0 to dynHigh-1
	    x = nz(filt[count])
	    y = -math.sin(2 * math.pi * count / dynHigh)
	    Sx := Sx + x
	    Sy := Sy + y
	    Sxx := Sxx + x * x
	    Sxy := Sxy + x * y
	    Syy := Syy + y * y
    
	imag = (dynHigh * Sxx - Sx*Sx > 0) and (dynHigh * Syy - Sy*Sy > 0) ? (dynHigh * Sxy - Sx*Sy) / math.sqrt((dynHigh * Sxx - Sx*Sx) * (dynHigh * Syy - Sy*Sy)) : 0
    
    //Compute the angle as an arctangent function and resolve ambiguity
    angle = real != 0 ? 90 - 180 * math.atan(imag/real) / math.pi : 0
    if (real < 0)
        angle := angle - 180
    
    //Do not allow the angle slope to go negative
    if (angle[1] - angle < 270 and angle < angle[1])
        angle := angle[1]
    
    //Limit DeltaAngle errors
    deltaAngle = angle - angle[1]
    if (deltaAngle <= 360/dynHigh)
        deltaAngle := deltaAngle[1]
    if (deltaAngle >= 360/dynLow)
        deltaAngle := deltaAngle[1]
    
    //Compute dominant cycle as summation of delta phases
    sumAngle = 0.0
    domCycle = 0
    for count = 0 to dynHigh
        sumAngle := sumAngle + deltaAngle[count]
        if sumAngle > 360
            domCycle := count
            break
    
    domCycle

// ----------
// Wrappers
// ----------

// @function Execute a particular Length Adaptation or Dominant Cycle Estimator from the list
// @param type Length Adaptation or Dominant Cycle Estimator type to use
// @param src Series to use
// @param len Reference lookback length
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @param chandeSDLen Lookback length of Standard deviation for Chande's Dynamic Length
// @param chandeSmooth Smoothing length of Standard deviation for Chande's Dynamic Length
// @param chandePower Exponent of the length adaptation for Chande's Dynamic Length (lower is smaller variation)
// @param preHP Use High Pass prefilter for the Estimators that support it (default)
// @param preSS Use Super Smoother prefilter for the Estimators that support it (default)
// @param preHP Use Hann Windowing prefilter for the Estimators that support it
// @returns Calculated period (float, not limited)
export doAdapt(string type, float src, int len, int dynLow, int dynHigh, int chandeSDLen = 5, int chandeSmooth = 10, float chandePower = 0.5, bool preHP = true, bool preSS = true, bool preHW = false) =>
    out = 0.0
    switch type
        "MESA MAMA Cycle" =>
            out := mamaPeriod(src, dynLow, dynHigh)
        "Pearson Autocorrelation" =>
            out := paPeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
        "DFT Cycle" =>
            out := dftPeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
        "Phase Accumulation" =>
            out := phasePeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
        =>
            out := la.doAdapt(type, src, len, dynLow, dynHigh, chandeSDLen, chandeSmooth, chandePower)
    
    out

// @function Execute a particular Dominant Cycle Estimator from the list
// @param type Dominant Cycle Estimator type to use
// @param src Series to use
// @param dynLow Lower bound for the dynamic length
// @param dynHigh Upper bound for the dynamic length
// @param preHP Use High Pass prefilter for the Estimators that support it (default)
// @param preSS Use Super Smoother prefilter for the Estimators that support it (default)
// @param preHP Use Hann Windowing prefilter for the Estimators that support it
// @returns Calculated period (float, not limited)
export doEstimate(string type, float src, int dynLow, int dynHigh, bool preHP = true, bool preSS = true, bool preHW = false) =>
    out = 0.0
    switch type
        "None" =>
            out := dynLow
        "MESA MAMA Cycle" =>
            out := mamaPeriod(src, dynLow, dynHigh)
        "Pearson Autocorrelation" =>
            out := paPeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
        "DFT Cycle" =>
            out := dftPeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
        "Phase Accumulation" =>
            out := phasePeriod(src, dynLow, dynHigh, preHP, preSS, preHW)
    
    out
