>> figure;, plot( HRSS(27310-270: 27310+329) )

> figure;, plot( HRSS(25450-270: 25450+329) )

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

0

100

200

300

400

500

600

>> figure;, plot( HRSS(22801-270: 22801+329) )

>> figure;, plot( HRSS(21410 -270: 21410 +329) )

>> [ploting_info] = MERLIN_DAMP_V1(HRSS(10000:1:end) ,64 ,256,1);

Note: this location

​

And this location

Both point to the same anomaly. One is pointing to the beginning, the other to the end

https://www.kaggle.com/inIT-OWL/high-storage-system-data-for-energy-optimization

The high storage system consists of 4 short conveyor belts (BLO, BHL, BHR, BRU) and 2 rails (HR, HL). The two conveyor belts in the middle (BHL, BHR) can be moved in the vertical by the rails, the other ones are fixed and they all have a size of 64cm x 8.5cm x 29.7cm.

​

Each conveyor belt has three induction sensors. The first one is 3.6cm from the left edge, the second one 26.6 cm from the left edge and the last sensor is 3.6cm from the right edge.

​

It uses a SPS with Codesys V3, which corresponds to IEC61131-Standard.

The high storage system transports one package between two spots, as you can see in https://www.youtube.com/watch?v=3o8PwyuwXXc&ab_channel=AlexandervB

The first run is the non-optimized run. The two conveyor belts in the middle are only moving vertical when they do not move the package horizontal. The second run is the optimized run. While the two conveyor belts in the middle are moving the package horizontal, they move vertical as well.

​

The generated data is split in four files. HRSSnormalstandard.csv contains normal runs without failures and not optimized.

​

HRSSnormaloptimized.csv containes optimized runs without failures.

HRSSanomalousstandard.csv contains runs with failures and not optimized.

And HRSSanomalousoptimized.csv contains optimized runs with failures.

​

The Label column in each file marks the rows with anomalies.

With these files you can test energy based optimization processes by using the normal non-optimized and normal optimized files.

�Furthermore you can test anomaly detection with the normal and anomaly files.

High Storage System Data for Energy Optimization

MADRID

Threshold (3 standard deviations from mean): 0.17589

Anomaly A: m = 64, loc = 51193, score = 0.48314

Anomaly B: m = 142, loc = 54229, score = 0.18922

Anomaly C: m = 252, loc = 58208, score = 0.18292

Dashed lines represent 3 different anomalies

Synthetic data

Madrid - Synthetic data

Train size = 50,000

TS length = 100,000

WarmStartDAMP

Elapsed time is 234.2799 seconds.

​

MADRID�Elapsed time is 264.7367 seconds.�MADRID semantically converged after 34.3180 seconds, at which point 1.56% of the total computation was completed.

​

After initialization (34.3180 seconds), MADRID has already converged, with all 193 predictions pointing to the correct anomalous positions. Among them, 58 point to Anomaly A, 89 point to Anomaly B, and 46 point to Anomaly C.

​

After WarmStartDAMP and MADRID completed their execution, out of 193 predictions, 83 point to anomaly A, 64 point to anomaly B, and 46 point to anomaly C.

Anomaly A

Anomaly B

Anomaly C

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

Madrid - Synthetic data

Train size = 50,000

TS length = 1 million

WarmStartDAMP

Elapsed time is 1958.8386 seconds.

​

MADRID�Elapsed time is 1988.9043 seconds.�MADRID semantically converged after 72.5268 seconds, at which point 1.55% of the total computation was completed.

​

After initialization (72.5268 seconds), MADRID has already converged, with all 193 predictions pointing to the correct anomalous positions. Among them, 58 point to Anomaly A, 89 point to Anomaly B, and 46 point to Anomaly C.

​

After WarmStartDAMP and MADRID completed their execution, out of 193 predictions, 83 point to anomaly A, 64 point to anomaly B, and 46 point to anomaly C.

Anomaly A

Anomaly B

Anomaly C

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

Comparison of MADRID and Warm-start DAMP

Parameter settings:

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

WarmStartDAMP

Elapsed time is 234.2799 seconds.

​

MADRID�Elapsed time is 264.7367 seconds.�MADRID semantically converged after 34.3180 seconds, at which point 1.56% of the total computation was completed.

Quality of solution

Percentage of processed pixels

TS length = 100,000

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

5

10

15

20

25

30

35

40

45

Warm-start DAMP

MADRID

semantically converged

Comparison of MADRID and Warm-start DAMP

Parameter settings:

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

WarmStartDAMP

Elapsed time is 234.2799 seconds.

​

MADRID�Elapsed time is 264.7367 seconds.�MADRID semantically converged after 34.3180 seconds, at which point 1.56% of the total computation was completed.

Quality of solution

Wall clock time

0

50

100

150

200

250

300

0

5

10

15

20

25

30

35

40

45

TS length = 100,000

Warm-start DAMP

MADRID

semantically converged

Comparison of MADRID and Warm-start DAMP

Parameter settings:

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

WarmStartDAMP

Elapsed time is 1958.8386 seconds.

​

MADRID�Elapsed time is 1988.9043 seconds.�MADRID semantically converged after 72.5268 seconds, at which point 1.55% of the total computation was completed.

Quality of solution

Percentage of processed pixels

Warm-start DAMP

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

MADRID

TS length = 1 million

semantically converged

Comparison of MADRID and Warm-start DAMP

Parameter settings:

minL = 64;

maxL = 256;

stepSize = 1;

trainTestSplit = 50000;

WarmStartDAMP

Elapsed time is 1958.8386 seconds.

​

MADRID�Elapsed time is 1988.9043 seconds.�MADRID semantically converged after 72.5268 seconds, at which point 1.55% of the total computation was completed.

Quality of solution

Wall clock time

TS length = 1 million

0

200

400

600

800

1000

1200

1400

1600

1800

2000

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Warm-start DAMP

MADRID

semantically converged

Telemanom

Train size = 50,000

L_s = 160

TS length = 1 million

Epoch 1/35

39864/39864 [==============================] - 482s 12ms/step - loss: 0.0277 - val_loss: 0.0095

Epoch 2/35

39864/39864 [==============================] - 536s 13ms/step - loss: 0.0142 - val_loss: 0.0080

Epoch 3/35

39864/39864 [==============================] - 525s 13ms/step - loss: 0.0131 - val_loss: 0.0081

Epoch 4/35

39864/39864 [==============================] - 577s 14ms/step - loss: 0.0123 - val_loss: 0.0078

Epoch 5/35

39864/39864 [==============================] - 594s 15ms/step - loss: 0.0120 - val_loss: 0.0078

Epoch 6/35

39864/39864 [==============================] - 597s 15ms/step - loss: 0.0120 - val_loss: 0.0084

Epoch 7/35

39864/39864 [==============================] - 505s 13ms/step - loss: 0.0118 - val_loss: 0.0077

Epoch 8/35

39864/39864 [==============================] - 572s 14ms/step - loss: 0.0117 - val_loss: 0.0077

Epoch 9/35

39864/39864 [==============================] - 532s 13ms/step - loss: 0.0115 - val_loss: 0.0084

Epoch 10/35

39864/39864 [==============================] - 560s 14ms/step - loss: 0.0115 - val_loss: 0.0077

Epoch 11/35

39864/39864 [==============================] - 550s 14ms/step - loss: 0.0116 - val_loss: 0.0080

Epoch 12/35

39864/39864 [==============================] - 593s 15ms/step - loss: 0.0114 - val_loss: 0.0078

Total time:

real 164m34.863s

user 92m53.027s

sys 233m33.119s

 start end

0 51120 51407

1 54150 54429

2 58210 58489

Anomaly found at row: 0

Anomaly found at row: 1

Anomaly found at row: 2

True positive: 3

False positive: 0

Top discord: 51263

Training time: 6653s = 110m53s

Testing time = 53m42s

Threshold (3 standard deviations from mean): 0.17589

Anomaly A: m = 64, loc = 51193, score = 0.48314

Anomaly B: m = 142, loc = 54229, score = 0.18922

Anomaly C: m = 252, loc = 58208, score = 0.18292

Dashed lines represent 3 different anomalies

top) The beginning and end of the test dataset. bottom) The three anomalies embedded into the synthetic data are visually obvious with human inspection.

minL = floor(1/3*DAMP_m)

maxL = 3*DAMP_m

202_UCR_Anomaly_CHARISfive_10411_10998_11028

DAMP_m = 32

​

minL = 10

maxL = 96

step = 1

success_num = 75

fail_num = 12

minL = floor(1/3*DAMP_m)

maxL = 3*DAMP_m

113_UCR_Anomaly_CIMIS44AirTemperature1_4000_5391_5392

DAMP_m = 24

​

minL = 8

maxL = 72

step = 1

DAMP_m = 24

success_num = 60

fail_num = 5

minL = floor(1/3*DAMP_m)

maxL = 3*DAMP_m

011_UCR_Anomaly_DISTORTEDECG1_10000_11800_12100

DAMP_m = 92

​

minL = 30

maxL = 276

step = 1

success_num = 194

fail_num = 53

DAMP_m = 92

 start end

0 6982 7238

1 55282 55351

2 66062 66131

3 103512 103861

Anomaly found at row: 0

Anomaly found at row: 3

True positive: 2

False positive: 2

Top discord: 7098

run_id = '2023-04-29_14.58.55’

train_test_split = 3000

l_s = 132

running time: 22m57.279s

Reported Anomalies:

 start end

0 6982 7238

1 55282 55351

2 66062 66131

3 103512 103861

Anomaly found at row: 0

Anomaly found at row: 3

True positive: 2

False positive: 2

Top discord: 7098

run_id = '2023-04-29_19.22.12'

train_test_split = 3000

l_s = 132

running time: 21m44.176s

Reported Anomalies:

start end

0 4382 4604

1 6982 7177

2 28612 28681

3 30642 30711

4 38552 38621

5 43802 43871

6 44012 44221

7 44292 44492

8 47232 47301

9 48352 48491

10 48632 48701

11 49542 49611

12 52412 52481

13 52622 52761

14 52972 53041

15 53252 53391

16 55282 55421

17 62492 62561

18 97842 97911

19 98402 98471

20 98682 98821

21 98962 99031

22 103582 103861

Anomaly found at row: 1

Anomaly found at row: 22

True positive: 2

False positive: 21

Top discord: 7028

run_id = '2023-04-27_19.59.33’

train_test_split = 3000

l_s = 132

running time: 23m46.267s

Telemanom

Reported Anomalies:

 start end

0 6944 7223

1 44044 44183

2 44324 44393

3 48384 48453

4 52654 52793

5 72534 72673

6 72814 72883

7 103544 103877

Anomaly found at row: 0

Anomaly found at row: 7

True positive: 2

False positive: 6

Top discord: 7098

run_id = '2023-04-29_22.42.05'

train_test_split = 3000

l_s = 164

running time: 26m38.567s

Reported Anomalies:

 start end

0 6944 7237

1 28574 28713

2 38584 38653

3 43764 43833

4 44044 44183

5 44324 44463

6 44604 44743

7 47194 47333

8 48314 48453

9 48664 48733

10 49504 49643

11 52374 52513

12 52654 52793

13 53004 53073

14 53284 53353

15 55244 55383

16 62454 62593

17 72534 72673

18 72814 72883

19 94304 94443

20 97804 97943

21 98364 98573

22 98644 98853

23 98994 99063

24 103544 103823

True positive: 2

False positive: 23

Top discord: 7098

run_id = '2023-04-29_21.31.59'

train_test_split = 3000

l_s = 164

running time: 36m43.806s

Reported Anomalies:

 start end

0 6524 6663

1 6944 7153

2 43764 43833

3 44044 44183

4 44324 44463

5 44604 44743

6 47194 47333

7 48384 48453

8 52374 52513

9 52654 52793

10 53004 53073

11 53284 53353

12 55244 55383

13 62454 62593

14 72534 72673

15 72814 72883

16 103544 103888

Anomaly found at row: 1

Anomaly found at row: 16

True positive: 2

False positive: 15

Top discord: 7032

run_id = '2023-04-29_21.31.59'

train_test_split = 3000

l_s = 164

running time: 38m59.005s

Reported Anomalies:

 start end

0 6950 7229

1 32500 32569

2 38590 38659

3 43770 43839

4 43980 44189

5 44260 44469

6 44610 44749

7 47200 47339

8 48320 48459

9 49510 49579

10 52380 52519

11 52660 52799

12 53010 53079

13 53290 53359

14 55250 55389

15 62460 62599

16 72540 72679

17 72820 72889

18 103550 103888

Anomaly found at row: 0

Anomaly found at row: 18

True positive: 2

False positive: 17

Top discord: 7098

run_id = '2023-04-29_20.27.33'

train_test_split = 3000

l_s = 100

running time: 19m40.819s

Reported Anomalies:

 start end

0 6950 7221

1 43770 43839

2 44050 44189

3 44330 44399

4 44610 44679

5 48390 48459

6 49510 49649

7 52380 52449

8 52660 52799

9 53010 53079

10 53290 53359

11 55250 55389

12 62460 62599

13 72540 72679

14 72820 72889

15 103550 103829

Anomaly found at row: 0

Anomaly found at row: 15

True positive: 2

False positive: 14

Top discord: 7085

run_id = '2023-04-29_20.08.45’

train_test_split = 3000

l_s = 100

running time: 18m36.445s

Reported Anomalies:

 start end

0 6950 7228

1 32500 32569

2 38590 38659

3 43770 43909

4 43980 44189

5 44330 44399

6 44610 44679

7 47200 47339

8 48320 48459

9 49510 49649

10 52380 52519

11 52660 52799

12 53010 53079

13 53290 53359

14 55250 55389

15 62460 62599

16 72540 72679

17 72820 72889

18 103550 103889

Anomaly found at row: 0

Anomaly found at row: 18

True positive: 2

False positive: 17

Top discord: 7085

run_id = '2023-04-29_19.49.31’

train_test_split = 3000

l_s = 100

running time: 17m58.921s

Reported Anomalies:

For reference only

function T= MultiScaleTSA

T =[];

Sr=0.062;

for i=1:9

if mod(i,2)

T = [ T sin(0:Sr:pi) ];

else

T = [ T sin(0:Sr:pi)/2 sin(0:Sr:pi)/2 ];

end

end

T = [ T sin(0:Sr:pi)/2 sin(0:Sr:pi)/2 sin(0:Sr:pi)/2 ];

for i=1:10

if mod(i,2)

T = [ T sin(0:Sr:pi) ];

else

T = [ T sin(0:Sr:pi)/2 sin(0:Sr:pi)/2 ];

end

end

end

Quality of solution

0

60,000

0%

100%

~15.5 hours

Pure brute force (NN)

Pure brute force (extrapolated)

~16 minutes

Pure brute force (NN)

(Zoom-in of the above)

Comparison of three brute-force algorithms

Test case: Sine wave of length 8192, with period 128 and anomaly position 4224-4352

Parameter settings:

minL = 128; % one cycle - one pi

maxL = 6*128; % six cycles

stepSize = 1;

train_test_split = 128*24;

Elapsed time for Warm-start DAMP is 1010.9067 seconds.

Elapsed time for DAMP is 1472.5665 seconds.

Elapsed time for Brute force (MASS) is 3913.77 seconds.

Elapsed time for Pure brute force (NN) is expected to be 56363.3 seconds.

Quality of solution

Wall clock time

Warm-start DAMP

Brute force (MASS)

Pure brute force (NN)

Comparison of three brute-force algorithms

Test case: Sine wave of length 8192, with period 128 and anomaly position 4224-4352

Parameter settings:

minL = 128; % one cycle - one pi

maxL = 6*128; % six cycles

stepSize = 1;

train_test_split = 128*24;

Elapsed time for Warm-start DAMP is 1010.9067 seconds.

Elapsed time for DAMP is 1472.5665 seconds.

Elapsed time for Brute force (MASS) is 3913.77 seconds.

Elapsed time for Pure brute force (NN) is expected to be 56363.3 seconds.

Quality of solution

Wall clock time

Warm-start DAMP

DAMP

Brute force (MASS)

Pure brute force (NN)

Comparison of three brute-force algorithms

Test case: Sine wave of length 8192, with period 128 and anomaly position 4224-4352

Parameter settings:

minL = 128; % one cycle - one pi

maxL = 6*128; % six cycles

stepSize = 1;

train_test_split = 128*24;

Elapsed time for Warm-start DAMP is 1010.9067 seconds.

Elapsed time for DAMP is 1472.5665 seconds.

Elapsed time for pure brute force is 3913.77 seconds.

0

500

1000

1500

2000

2500

3000

3500

4000

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Quality of solution

Wall clock time

Warm-start DAMP

DAMP

pure brute force

function testAnomaly

T=2*rem(1:1000,2)/2;

T([1:50:1000])= 5;

T(430:431)=[];

[matrixProfile, profileIndex, motifIndex, discordIndex] = interactiveMatrixProfile( T,49 ); %48 no, 49 yes

max(matrixProfile)

figure;

T=2*rem(1:1000,2)/2;

for i = 1:100:1000

T(i:i+50)= T(i:i+50)+100;

end

T(453:499) = randn(size(T(453:499)))/5;

[matrixProfile, profileIndex, motifIndex, discordIndex] =interactiveMatrixProfile( T,50 ); % 49 yes, 50 no

max(matrixProfile)

end

​

>> figure;, MERLIN_DAMP_V1(melbourne7(66820-(24*7*30):66820 +(7 * 6 * 24),1) ,8,1*24,1)

As shown in Figure X.top we have an eight-minute long dataset containing lead III of an ECG  (footnote) [xx]. This is the only data that MADRID is allowed to see for this experiment. However, cardiologist Dr. <blinded> is also able to examine other telemetry that was recorded in parallel and the original clinical notes and was thus able to produce the definitive ground truth. The ground truth is that this dataset has two bradycardia arrhythmias, one that is visually obvious and one that is very subtle, perhaps undetectable to most cardiologists that did not also have access to the out-of-band data.

Can MADRID find these anomalies? Figure X shows the answer. What is remarkable is that the more subtle anomaly only shows up at a fairly limited range of lengths, from 49 to 63.

This clearly motivates the use of MADRID. Even though both anomalies have the same cardiological cause, they manifested themselves at different scales, and these scales that would be impossible to know ahead of time.

 (footnote) The data in fig X.top has a curious property that all the peaks have exactly the same value. This was due to a mechanical limitation of the sensor used, and is inconsequential for the task at hand

[xx] Moody GB, Mark RG. A Database to Support Development and Evaluation of Intelligent Intensive Care Monitoring. Computers in Cardiology 23:657–660 (1996).

​

Figure X top) An eight-minute-long segment of an ECG. The first eight seconds (about nine beats) are used as training data. center) MADRID discovered every anomaly from length 32 to 128, which are clustered to into just two locations. bottom) The two discovered anomalies are both bradycardia arrhythmias.

450m.mat

Supporting data for Bradycardia arrhythmia example

Supporting data for Bradycardia arrhythmia example

Supporting data for Bradycardia arrhythmia example

Supporting data for Bradycardia arrhythmia example

ploting_info =

​

1.0e+04 *

​

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0.0032 1.0382 0.0000

0.0033 1.0382 0.0000

0.0034 1.0387 0.0000

0.0035 1.0386 0.0000

0.0036 1.0387 0.0000

0.0037 1.0386 0.0000

0.0038 1.0385 0.0000

0.0039 1.0384 0.0000

0.0040 1.0383 0.0000

0.0041 1.0382 0.0000

0.0042 1.0381 0.0000

0.0043 1.0380 0.0000

0.0044 1.0379 0.0000

0.0045 1.0378 0.0000

0.0046 1.0377 0.0000

0.0047 1.0376 0.0000

0.0048 1.0375 0.0000

0.0049 2.5390 0.0000

0.0050 2.5389 0.0000

0.0051 2.5388 0.0000

0.0052 2.5387 0.0000

0.0053 2.5386 0.0000

0.0054 2.5385 0.0000

0.0055 2.5384 0.0000

0.0056 2.5383 0.0000

0.0057 2.5382 0.0000

0.0058 2.5381 0.0000

0.0059 2.5378 0.0000

0.0060 2.5377 0.0000

0.0061 2.5377 0.0000

0.0062 2.5376 0.0000

0.0063 2.5375 0.0000

0.0064 1.0382 0.0000

0.0065 1.0381 0.0000

0.0066 1.0380 0.0000

0.0067 1.0379 0.0000

0.0068 1.0378 0.0000

0.0069 1.0377 0.0000

0.0070 1.0376 0.0000

0.0071 1.0375 0.0000

0.0072 1.0374 0.0000

0.0073 1.0373 0.0000

0.0074 1.0372 0.0000

0.0075 1.0371 0.0000

0.0076 1.0370 0.0000

0.0077 1.0369 0.0000

0.0078 1.0368 0.0000

0.0079 1.0367 0.0000

0.0080 1.0366 0.0000

0.0081 1.0365 0.0000

0.0082 1.0364 0.0000

0.0083 1.0363 0.0000

0.0084 1.0362 0.0000

0.0085 1.0361 0.0000

0.0086 1.0360 0.0000

0.0087 1.0359 0.0000

0.0088 1.0358 0.0000

0.0089 1.0357 0.0000

0.0090 1.0356 0.0000

0.0091 1.0355 0.0000

0.0092 1.0354 0.0000

0.0093 1.0353 0.0000

0.0094 1.0352 0.0000

0.0095 1.0351 0.0000

0.0096 1.0350 0.0000

0.0097 1.0349 0.0000

0.0098 1.0348 0.0000

0.0099 1.0347 0.0000

0.0100 1.0346 0.0000

0.0101 1.0345 0.0000

0.0102 1.0344 0.0000

0.0103 1.0342 0.0000

0.0104 1.0341 0.0000

0.0105 1.0339 0.0000

0.0106 1.0338 0.0000

0.0107 1.0335 0.0000

0.0108 1.0288 0.0000

0.0109 1.0287 0.0000

0.0110 1.0312 0.0000

0.0111 1.0311 0.0000

0.0112 1.0312 0.0000

0.0113 1.0311 0.0000

0.0114 1.0333 0.0000

0.0115 1.0332 0.0000

0.0116 1.0331 0.0001

0.0117 1.0330 0.0001

0.0118 1.0330 0.0001

0.0119 1.0329 0.0001

0.0120 1.0328 0.0001

0.0121 1.0327 0.0001

0.0122 1.0326 0.0001

0.0123 1.0325 0.0001

0.0124 1.0324 0.0001

0.0125 1.0323 0.0001

0.0126 1.0322 0.0001

0.0127 1.0321 0.0001

0.0128 1.0320 0.0001

This is output used to make the plots

​

​

Column 1 is subsequence length. Lengths 1 to 31 are all zeros, they are just place keepers

​

​

Column 2 is the location of the discord

​

​

Column 3 is the value of the discord, in length normalized form

​

​

Supporting data for Bradycardia arrhythmia example

function [ploting_info] = MERLIN_DAMP_V1(T,minL,maxL,step)

k = 1; % how many anomalies do you want (at every length;

train_length = 1000; % How long of a prefix do you want to consider as training data

if ~CheckInput(T,minL,maxL), return;, end % Make sure parameters are sensible

MERLIN_STEM_PLOTindex =[]; % length of subsequence used

MERLIN_STEM_PLOTdiscordLoc =[];

MERLIN_STEM_PLOTdiscordVal =[];

for i = minL:step:maxL

disp(['Step ', num2str(i), ' of ',num2str(maxL),]);

MERLIN_STEM_PLOTindex(i) = i; % For a subsequence of length i

resutls = DAMP_topK(T,train_length,i,k,false);

MERLIN_STEM_PLOTdiscordLoc(i,:)= resutls(:,2)';

MERLIN_STEM_PLOTdiscordVal(i,:)= resutls(:,1)' * (1/(2*sqrt(i)));

end

hold on;

for i = 1 : k % loop over the top k anomalies

stem3( MERLIN_STEM_PLOTindex(:), MERLIN_STEM_PLOTdiscordLoc(:,i) , MERLIN_STEM_PLOTdiscordVal(:,i) )

end

set(gca,'Ylim',[1,length(T)])

set(gca,'Xlim',[minL,maxL])

set(gca,'xdir', 'reverse');

T_normalized_for_plot = (T-min(T))/max(T);

T_normalized_for_plot = T_normalized_for_plot * max(MERLIN_STEM_PLOTdiscordVal(:,1));

plot3( ones(size([1:length(T)]))* maxL , [1:length(T)],T_normalized_for_plot,'b' );

plot3( ones(size([1:train_length]))* maxL , [1:train_length],T_normalized_for_plot(1:train_length),'r' );

view([59 68]);

[ploting_info] = [MERLIN_STEM_PLOTindex', MERLIN_STEM_PLOTdiscordLoc , MERLIN_STEM_PLOTdiscordVal];

end

function SensibleParameters = CheckInput(T,minL,maxL)

SensibleParameters = false; % Assume the parameters are bad, then test them

if (maxL * 10) > length(T)

disp('time series must be at least 10 times the length of maxL ');

elseif minL < 3

disp('minL must be at least 3 ');

elseif (maxL - minL) < 0

disp('minL must be shorter than maxL ');

else

SensibleParameters = true; % paramters passed all our tests

end

end

function Results = DAMP(T,CurrentIndex,SubsequenceLength,discord_num,enable_output)

% This is a special Matrix Profile, it only looks left (backwards in

% time)

Left_MP = zeros(size(T));

% Initialization

% The best discord score so far

best_so_far = -inf;

% A Boolean vector where 1 means execute the current iteration and 0

% means skip the current iteration

bool_vec = ones(1,length(T));

% Lookahead indicates how long the algorithm has a delay

lookahead = 2^nextpow2(16*SubsequenceLength);

% Handle the prefix to get a relatively high best so far discord score

% Prefix for top k

for i = CurrentIndex : CurrentIndex+16*SubsequenceLength

% Use the brute force for the left Matrix Profile value

if i+SubsequenceLength-1 > length(T)

break

end

Left_MP(i) = min( real(MASS_V2(T(1:i), T(i:i+SubsequenceLength-1))));

end

Left_MP_copy = Left_MP;

for k= 1:discord_num

[best_so_far,idx_max]=max(Left_MP_copy);

discord_start = max(1,idx_max - floor(SubsequenceLength*0.5));

discord_end = max(1+floor(SubsequenceLength*0.5),idx_max + floor(SubsequenceLength*0.5));

Left_MP_copy(discord_start:discord_end)=-inf;

end

% Remaining test data except for the prefix

for i = (CurrentIndex+(16*SubsequenceLength)+1) : (length(T)-SubsequenceLength+1)

% Skip the current iteration if the corresponding boolean value is

% 0, otherwise execute the current iteration

if ~bool_vec(i)

% We subtract a very small number here to avoid the pruned

% subsequence having the same discord score as the real discord

Left_MP(i) = Left_MP(i-1)-0.00001;

continue

end

% Use the brute force for the left Matrix Profile value

if i+SubsequenceLength-1 > length(T)

break

end

% Initialization for classic DAMP

% Approximate leftMP value for the current subsequence

approximate_distance = inf;

% X indicates how long a time series to look backwards

X = 2^nextpow2(8*SubsequenceLength);

% flag indicates if it is the first iteration of DAMP

flag = 1;

% expansion_num indicates how many times the search has been

% expanded backward

expansion_num = 0;

query = T(i:i+SubsequenceLength-1);

% Classic DAMP

while approximate_distance >= best_so_far

% Case 1: Execute the algorithm on the time series segment

% farthest from the current subsequence

% Arrived at the beginning of the time series

if i-X+1+(expansion_num * SubsequenceLength) < 1

approximate_distance = min( real(MASS_V2(T(1:i), query)));

Left_MP(i) = approximate_distance;

% Update the best discord so far

if approximate_distance > best_so_far

% The current subsequence is the best discord so far

best_so_far = approximate_distance;

Left_MP_copy = Left_MP;

for k= 1:discord_num

[best_so_far,idx_max]=max(Left_MP_copy);

discord_start = max(1,idx_max - floor(SubsequenceLength*0.5));

discord_end = max(1+floor(SubsequenceLength*0.5),idx_max + floor(SubsequenceLength*0.5));

Left_MP_copy(discord_start:discord_end)=-inf;

end

end

break

else

if flag == 1

% Case 2: Execute the algorithm on the time series

% segment closest to the current subsequence

flag = 0;

approximate_distance = min( real(MASS_V2(T(i-X+1:i), query)));

else

% Case 3: All other cases

X_start = i-X+1+(expansion_num * SubsequenceLength);

X_end = i-(X/2)+(expansion_num * SubsequenceLength);

approximate_distance = min( real(MASS_V2(T(X_start:X_end), query)));

end

if approximate_distance < best_so_far

% If a value less than the current best discord score

% exists on the distance profile, stop searching

Left_MP(i) = approximate_distance;

break

else

% Otherwise expand the search

X = 2*X;

expansion_num = expansion_num+1;

end

end % end if

end % end while

% If lookahead is 0, then it is a pure online algorithm with no

% pruning

if lookahead ~= 0

% Perform forward MASS for pruning

% The index at the beginning of the forward mass should be

% avoided in the exclusion zone

start_of_mass = i+SubsequenceLength;

if start_of_mass > length(T)

start_of_mass = length(T);

end

end_of_mass = start_of_mass + lookahead - 1;

if end_of_mass > length(T)

end_of_mass = length(T);

end

% The length of lookahead should be longer than that of the

% query

if (end_of_mass - start_of_mass + 1) > SubsequenceLength

distance_profile = real(MASS_V2(T(start_of_mass:end_of_mass), query));

% Find the subsequence indices less than the best so far

% discord score

dp_index_less_than_BSF = find((distance_profile<best_so_far)==1);

% Converting indexes on distance profile to indexes on time

% series

ts_index_less_than_BSF = dp_index_less_than_BSF + start_of_mass - 1;

% update the Boolean vector

bool_vec(ts_index_less_than_BSF) = 0;

end

end

end % end for

Results = zeros(discord_num,2);

% Print pruning rate

PV = bool_vec(CurrentIndex:(length(T)-SubsequenceLength+1));

PR = (length(PV)-sum(PV))/(length(PV));

if enable_output

disp(strcat("Pruning Rate: ",num2str(PR)));

end

% Print top K discords

Left_MP_copy = Left_MP;

for k= 1:discord_num

[val,loc]=max(Left_MP_copy);

if val == 0

if enable_output

disp(strcat("Only ",num2str(k-1)," discords are found."));

end

break

end

if enable_output

disp(strcat("Predicted discord score/position","(top ",num2str(k),"): ",num2str(val),"/",num2str(loc)));

end

Results(k,:)=[val,loc];

discord_start = max(1,loc - floor(SubsequenceLength*0.5));

discord_end = max(1+ floor(SubsequenceLength*0.5),loc + floor(SubsequenceLength*0.5));

Left_MP_copy(discord_start:discord_end)=-inf;

end

if enable_output

% Create the plot

figure; hold on; plot(Left_MP,'b'); plot(zscore(T)-2,'r');

end

end

function [dist] = MASS_V2(x, y)

%x is the data, y is the query

m = length(y);

n = length(x);

%compute y stats -- O(n)

meany = mean(y);

sigmay = std(y,1);

%compute x stats -- O(n)

meanx = movmean(x,[m-1 0]);

sigmax = movstd(x,[m-1 0],1);

y = y(end:-1:1);%Reverse the query

y(m+1:n) = 0; %aappend zeros

%The main trick of getting dot products in O(n log n) time

X = fft(x);

Y = fft(y);

Z = X.*Y;

z = ifft(Z);

dist = 2*(m-(z(m:n)-m*meanx(m:n)*meany)./(sigmax(m:n)*sigmay));

dist = sqrt(dist);

end

​

This is the exact code I used

This is the exact call I made

​

>> [ploting_info] = MERLIN_DAMP_V1(val(1,1:1:1:60000) ,32,128,1);

Supporting data for Bradycardia arrhythmia example