stelmo / linearsegmentation.jl Goto Github PK

# External
using LinearSegmentation;
using StableRNGs;
using UnicodePlots;

## Constants & Configuration

oRng = StableRNG(123);

## Functions

function Conv1D( vA :: Vector{T}, vB :: Vector{T}; convMode :: String = "full" ) :: Vector{T} where {T <: Real}

    lenA = length(vA);
    lenB = length(vB);

    if (convMode == "full")
        startIdx    = 1;
        endIdx      = lenA + lenB - 1;
    elseif (convMode == "same")
        startIdx    = 1 + floor(Int, lenB / 2);
        endIdx      = startIdx + lenA - 1;
    elseif (convMode == "valid")
        startIdx    = lenB;
        endIdx      = lenA;
    end

    vO = zeros(T, lenA + lenB - 1);

    for idxB in 1:lenB
        @simd for idxA in 1:lenA
            @inbounds vO[idxA + idxB - 1] += vA[idxA] * vB[idxB];
        end
    end

    return vO[startIdx:endIdx];
end

## Parameters
# Data
numSamples    = 500;
ampFiltSize   = 20;
phaseFiltSize = 50;

vSeg = [1, 151, 301, 401, 501];

# Model
minSegLen = 30.0;
maxRmse   = 0.12;

## Load / Generate Data

vAmp    = rand(oRng, numSamples);
vAmp    = 0.2 * Conv1D(vAmp, ones(ampFiltSize) / ampFiltSize; convMode = "same");
vPhase  = 0.2 * rand(oRng, numSamples);
vPhase  = Conv1D(vPhase, ones(phaseFiltSize) / phaseFiltSize; convMode = "same");
vPhase  = cumsum(vPhase);

vX = LinRange(0, numSamples - 1, numSamples);

vC = vAmp .* cos.(2 * pi * vPhase);
vL = zeros(numSamples);
vL[vSeg[1]:(vSeg[2] - 1)] .= 0;
vL[vSeg[2]:(vSeg[3] - 1)] .= 1;
vL[vSeg[3]:(vSeg[4] - 1)] .= collect(LinRange(0.5, 1.0, 100));
vL[vSeg[4]:(vSeg[5] - 1)] .= collect(LinRange(1.0, 0.4, 100));

vY = vC .+ vL;

## Display Data

ii = 1;
vIdx = vSeg[ii]:(vSeg[ii + 1] - 1)

hP = scatterplot(vX[vIdx], vY[vIdx], width = 90, height = 8, xlim = (vX[1], vX[end]), ylim = (minimum(vY), maximum(vY)));

for ii in 2:(length(vSeg) - 1)
    local vIdx = vSeg[ii]:(vSeg[ii + 1] - 1)
    scatterplot!(hP, vX[vIdx], vY[vIdx]);
end

title!(hP, "Input Data");
xlabel!(hP, "Index");
ylabel!(hP, "Value");
display(hP);

## Analysis
segs = shortest_path(vX, vY; min_segment_length = minSegLen, fit_function = :rmse, fit_threshold = maxRmse);

# Remove the 1st item which is shared
for ii in 2:length(segs)
    deleteat!(segs[ii][1], 1);
end

## Display Results
hP = scatterplot(segs[1][1], vY[segs[1][1]], width = 90, height = 8, xlim = (vX[1], vX[end]), ylim = (minimum(vY), maximum(vY)));

for ii in 2:length(segs)
    scatterplot!(hP, segs[ii][1], vY[segs[ii][1]]);
end

title!(hP, "Linear Segmentation");
xlabel!(hP, "Index");
ylabel!(hP, "Value");
display(hP);