11---
22jupyter :
33 jupytext :
4+ formats : ipynb,md
45 text_representation :
56 extension : .md
67 format_name : markdown
78 format_version : ' 1.3'
8- jupytext_version : 1.14.4
9+ jupytext_version : 1.14.5
910 kernelspec :
10- display_name : env_qolmat
11+ display_name : Python 3 (ipykernel)
1112 language : python
12- name : env_qolmat
13+ name : python3
1314---
1415
1516``` python
@@ -28,33 +29,24 @@ from math import pi
2829from qolmat.utils import plot, data
2930from qolmat.imputations.rpca.rpca_pcp import RPCAPCP
3031from qolmat.imputations.rpca.rpca_noisy import RPCANoisy
32+ from qolmat.utils.data import generate_artificial_ts
3133```
3234
3335** Generate synthetic data**
3436
3537``` python
3638n_samples = 1000
39+ periods = [100 , 20 ]
40+ amp_anomalies = 0.5
41+ ratio_anomalies = 0.05
42+ amp_noise = 0.1
3743
38- mesh = np.arange(n_samples)
39- X_true = np.zeros(n_samples)
40- A_true = np.zeros(n_samples)
41- E_true = np.zeros(n_samples)
42- p1 = 100
43- p2 = 20
44- X_true = 1 + np.sin(2 * pi * mesh / p1) + np.sin(2 * pi * mesh / p2)
45- noise = np.random.uniform(size = n_samples)
46- amplitude_A = .5
47- freq_A = .05
48- A_true = amplitude_A * np.where(noise < freq_A, - np.log(noise), 0 ) * (2 * (np.random.uniform(size = n_samples) > .5 ) - 1 )
49- amplitude_E = .1
50- E_true = amplitude_E * np.random.normal(size = n_samples)
51-
52- signal = X_true + E_true
53- signal[A_true != 0 ] = A_true[A_true != 0 ]
54- signal = signal.reshape(- 1 , 1 )
44+ X_true, A_true, E_true = generate_artificial_ts(n_samples, periods, amp_anomalies, ratio_anomalies, amp_noise)
45+
46+ signal = X_true + A_true + E_true
5547
5648# Adding missing data
57- signal[5 :20 , 0 ] = np.nan
49+ signal[5 :20 ] = np.nan
5850```
5951
6052``` python
@@ -73,7 +65,7 @@ plt.plot(E_true)
7365
7466ax = fig.add_subplot(4 , 1 , 4 )
7567ax.title.set_text(" Corrupted signal" )
76- plt.plot(signal[:, 0 ] )
68+ plt.plot(signal)
7769
7870plt.show()
7971```
@@ -82,26 +74,78 @@ plt.show()
8274
8375``` python
8476%% time
85-
8677rpca_pcp = RPCAPCP(period = 100 , max_iter = 5 , mu = .5 , lam = 1 )
87- X = rpca_pcp.fit_transform(signal)
88- corruptions = signal - X
78+ X, A = rpca_pcp.decompose_rpca_signal(signal)
79+ imputed = signal - A
80+ ```
81+
82+ ``` python
83+ fig = plt.figure(figsize = (12 , 4 ))
84+ plt.plot(X, color = " black" )
85+ plt.plot(imputed)
8986```
9087
9188## Temporal RPCA
9289
9390``` python
91+ %% time
9492rpca_noisy = RPCANoisy(period = 10 , tau = 2 , lam = 0.3 , list_periods = [10 ], list_etas = [0.01 ], norm = " L2" )
95- X = rpca_noisy.fit_transform(signal)
96- corruptions = signal - X
97- plot.plot_signal([signal[:,0 ], X[:,0 ], corruptions[:, 0 ]])
93+ X, A = rpca_noisy.decompose_rpca_signal(signal)
94+ ```
95+
96+ ``` python
97+ fig = plt.figure(figsize = (12 , 4 ))
98+ plt.plot(X, color = " black" )
99+ plt.plot(imputed)
100+ ```
101+
102+ ``` python
103+
104+ ```
105+
106+ ``` python
107+ %% time
108+ signal_toy = np.array([[1 , 2 ], [np.nan, np.nan]])
109+ rpca_noisy = RPCANoisy(tau = 0 , lam = 1 , norm = " L2" , do_report = True )
110+ X, A = rpca_noisy.decompose_rpca_signal(signal_toy)
111+ ```
112+
113+ ``` python
114+ print (X)
115+ print (A)
116+ ```
117+
118+ ``` python
119+ %% time
120+ signal_toy = np.array([[1 , 2 ], [np.nan, np.nan]])
121+ rpca_pcp = RPCAPCP(lam = 1e3 )
122+ X, A = rpca_pcp.decompose_rpca_signal(signal_toy)
123+ ```
124+
125+ ``` python
126+ X
127+ ```
128+
129+ ``` python
130+ A
131+ ```
132+
133+ ``` python
134+ np.log(10 ) / np.log(1.1 )
135+ ```
136+
137+ ``` python
138+ X = np.array([[1 , 2 ], [4 , 4 ], [4 , 3 ]])
139+ # Omega = np.array([[True, False], [True, True], [False, True]])
140+ Omega = np.array([[True , True ], [True , True ], [True , True ]])
141+ rpca_noisy = RPCANoisy(period = 2 , max_iter = 200 , tau = .5 , lam = 1 , do_report = True )
142+ M_result, A_result, U_result, V_result = rpca_noisy.decompose_rpca_L2(
143+ X, Omega = Omega, lam = 1 , tau = .5 , rank = 2
144+ )
98145```
99146
100147``` python
101- rpca_noisy = RPCANoisy(period = 10 , tau = 2 , lam = 0.3 , list_periods = [], list_etas = [], norm = " L2" )
102- X = rpca_noisy.fit_transform(signal)
103- corruptions = signal - X
104- plot.plot_signal([signal[:,0 ], X[:,0 ], corruptions[:, 0 ]])
148+ M_result
105149```
106150
107151``` python
0 commit comments