projects / epclust.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
'update'
[epclust.git] / temp / WerDist.R
CommitLineData
b7cd987d
BA		 1#Ligne 83 dans epclust/R/computeWerDists.R :
2# TODO:
3# [X,period,scale,coix] =
4#wavelet(x(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);%#ok
5#[Y,period,scale,coiy] = wavelet(y(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
6#
7#%Smooth X and Y before truncating! (minimize coi)
8#%sinv=1./(scale');
9#%
10#%
11#%sX=smoothwavelet(sinv(:,ones(1,nx)).*(abs(X).^2),dt,period,Args.Dj,scale);
12#%sY=smoothwavelet(sinv(:,ones(1,ny)).*(abs(Y).^2),dt,period,Args.Dj,scale);
13#%
14#%Wxy=X.*conj(Y);
15#%
16#%% ----------------------- Wavelet coherence ---------------------------------
17#%sWxy=smoothwavelet(sinv(:,ones(1,n)).*Wxy,dt,period,Args.Dj,scale);
18#%Rsq=abs(sWxy).^2./(sX.*sY);
19#
20
21
22#https://github.com/grinsted/wavelet-coherence/blob/master/wtc.m
23
24
25library(Rwave) # CWT
26library(cluster) # pam
27#library(flexclust) # kcca
28
29
30 toCWT <- function(X, sw= 0, tw= 0, swabs= 0,
31 nvoice= 12, noctave= 5,
32 s0= 2, w0= 2*pi, lt= 24, dt= 0.5,
33 spectra = FALSE, smooth = TRUE,
34 scaled = FALSE,
35 scalevector)
36 { noctave <- adjust.noctave(lt, dt, s0, tw, noctave)
37 if(missing(scalevector))
38 scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
39 res <- lapply(1:nrow(X), function(n)
40 { tsX <- ts( X[n,] )
41 tsCent <- tsX - mean(tsX)
42 if(scaled) tsCent <- ts(scale(tsCent))
43 tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0)
44 tsCent.cwt
45 } )
46 if( spectra ) res <- lapply(res, function(l) Mod(l)^2 )
47 if( smooth ) res <- lapply(res, smCWT, swabs = swabs,
48 tw = tw, dt = dt,
49 scalevector = scalevector)
50 resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]),
51 length(res)))
52 for( l in 1:length(res) ) resArray[ , , l] <- res[[l]]
53 resArray
54 }
55
56
57 # ===============================================================
58
59 smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0,
60 nvoice= 12, noctave= 2, s0= 2, w0= 2*pi,
61 lt= 24, dt= 0.5, scalevector )
62 {
63# noctave <- adjust.noctave(lt, dt, s0, tw, noctave)
64# scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
65 wsp <- Mod(CWT)
66 smwsp <- smooth.matrix(wsp, swabs)
67 smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
68 smsmwsp
69 }
70
71
72 # ===============================================================
73
74 toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm")
75{ x2 <- spline(x, n = 2^ceiling( log(length(x), 2) ),
76 method = 'natural')$y
77 Dx2 <- wd(x2, family = family, filter.number = filter.number)$D
78 Dx2
79}
80
81 # ===============================================================
82
83 contrib <- function(x)
84 { J <- log( length(x)+1, 2)
85 nrj <- numeric(J)
86 t0 <- 1
87 t1 <- 0
88 for( j in 1:J ) {
89 t1 <- t1 + 2^(J-j)
90 nrj[j] <- sqrt( sum( x[t0:t1]^2 ) )
91 t0 <- t1 + 1
92 }
93 return(nrj)
94 }
95
96
97 # ========================================= distance for coh ===
98
99 coherence <- function( x, y)
100 { J <- log(length(x) + 1, 2)
101 t0 <- 1
102 sg2_x <- 0
103 sg2_y <- 0
104 sg_xy <- 0
105 for(j in 0:(J - 1))
106 { t1 <- t0 + 2^(J - j)/2 - 1
107 tt <- t0:t1
108 sg2_x <- sg2_x + mean(x[t0:t1]^2)
109 sg2_y <- sg2_y + mean(y[t0:t1]^2)
110 sg_xy <- sg_xy + mean(x[t0:t1] * y[t0:t1])
111 t0 <- t1 + 1
112 }
113 res <- sg_xy^2 / sg2_x / sg2_y
114 res
115 }
116
117
118 vect2mat <- function(vect){
119 vect <- as.vector(vect)
120 matrix(vect[-(1:2)], delta, lscvect)
121 }
122
123
124 # ========================================= # myimg for graphics
125 jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF",
126 "cyan", "#7FFF7F", "yellow",
127 "#FF7F00", "red", "#7F0000"))
128
129 myimg <- function(MAT, x = 1:nrow(MAT), y = 1:col(MAT), ... )
130 filled.contour( x = x, y = y, z = MAT,
131 xlab= 'Time', ylab= 'scale',
132 color.palette = jet.colors,
133 ... )
134
135
136#TODO: [plus tard] alternative à sowa (package disparu) : cwt..
137#source("sowas-superseded.r") # auxiliary CWT functions
138
139
140cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi){
141
142 if (class(ts)!="ts"){
143
144 cat("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n")
145
146 }
147 else{
148
149 t=time(ts)
150 dt=t[2]-t[1]
151
152 s0unit=s0/dt*w0/(2*pi)
153 s0log=as.integer((log2(s0unit)-1)*nvoice+1.5)
154
155 if (s0log<1){
156 cat(paste("# s0unit = ",s0unit,"\n",sep=""))
157 cat(paste("# s0log = ",s0log,"\n",sep=""))
158 cat("# s0 too small for w0! \n")
159 }
160 totnoct=noctave+as.integer(s0log/nvoice)+1
161
162 totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0)
163
164 ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)]
165
166 #Normalization
167 sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
168 smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE)
169
170 ts.cwt*smat
171
172 }
173
174}
175
176#####################################
177# WSP
178#####################################
179
180wsp <- function(ts,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,siglevel=0.95,critval=NULL,nreal=1000,arealsiglevel=0.9,kernel=0,markt=-999,marks=-999,logscale=FALSE,plot=TRUE,units="",device="screen",file="wsp",color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=3){
181
182 if (class(ts)!="ts"){
183
184 cat("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n")
185
186 }
187 else{
188
189 if (sw!=0 & swabs==0)
190 swabs <- as.integer(sw*nvoice)
191 if (swabs!=0 & sw==0)
192 sw <- swabs/nvoice
193
194 sllist <- checkarealsiglevel(sw,tw,w0,arealsiglevel,siglevel,0)
195 arealsiglevel <- sllist$arealsiglevel
196 siglevel <- sllist$siglevel
197
198 at <- NULL
199
200 t <- time(ts)
201 dt <- deltat(ts)
202 s0rem <- s0
203 s0 <- adjust.s0(s0,dt)
204 dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1
205
206 noctave <- adjust.noctave(length(ts),dt,s0,tw,noctave)
207 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
208 tsanom <- ts-mean(ts)
209
210 #WAVELET TRANSFORMATION
211 ts.cwt <- cwt.ts(tsanom,s0,noctave,nvoice,w0)
212
213 #SMOOTHING
214 wsp <- smooth.matrix(Mod(ts.cwt)^2,swabs)
215 smwsp <- smooth.time(wsp,tw,dt,scalevector)
216
217 #POINTWISE SIGNIFICANCE TEST
218 if (is.null(critval)==FALSE){ # is critval empty?
219 if (dim(critval)[2]!=dim(smwsp)[2]){ # is critval of the wrong dimension?
220 if (siglevel[1]!=0 & nreal!=0) critval <-
221 criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal)
222 #critval is of wrong dimension and siglevel and nreal are given
223 else {
224 critval <- NULL # no test possible
225 arealsiglevel <- 0
226 cat("# dimension of critval is wrong \n")
227 cat("# areawise test only possible with pointwise test \n")
228 }
229 }
230 }
231 else{ # critval is empty, but nreal or siglevel is given
232 if (siglevel[1]!=0 & nreal!=0) critval <-
233 criticalvaluesWSP(tsanom,s0,noctave,nvoice,w0,swabs,tw,siglevel,nreal)
234 else {
235 critval <- NULL
236 arealsiglevel <- 0
237 cat("# areawise test only possible with pointwise test \n")
238 }
239 }
240
241 #AREAL SIGNIFICANCE TEST
242 if (arealsiglevel!=0){
243 v <- critval[1,]
244 v[v==0] <- 9999
245 cvmat <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE)
246 atest <- arealtest(smwsp/cvmat,dt,s0,noctave,nvoice,w0,swabs,tw,siglevel,arealsiglevel,kernel,0)
247 at <- atest$at
248 kernel <- atest$kernel
249 }
250
251 if (s0rem<s0){
252 smwsp <- addvalues(nvoice,dnoctave,smwsp,NA)
253 critval <- addvalues(nvoice,dnoctave,critval,1)
254 #at <- addvalues(nvoice,dnoctave,at,NA)
255 noctave <- noctave+dnoctave
256 s0 <- s0/2^dnoctave
257 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
258 }
259
260 #PARAMETERS
261 wclist <-
262 list(modulus=smwsp,phase=NULL,time=t,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,scales=scalevector,critval=critval,at=at,kernel=kernel)
263
264 class(wclist) <- "wt"
265
266 #PLOTTING
267 if (plot)
268 plot(wclist,markt,marks,NULL,NULL,logscale,FALSE,units,"Wavelet Power Spectrum",device,file,FALSE,color,pwidth,pheight,labsc,labtext,sigplot)
269
270 wclist
271
272 }
273
274}
275
276#####################################
277# WCO
278#####################################
279
280wco <-
281 function(ts1,ts2,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,siglevel=0.95,arealsiglevel=0.9,kernel=0,markt=-999,marks=-999,sqr=FALSE,phase=TRUE,plot=TRUE,units="",device="screen",file="wcoh",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=3){
282
283 if (class(ts1)!="ts"){
284
285 cat("# This function needs two time series objects as input. You may construct them by using the function ts(data,start,deltat). Try '?ts' for help.\n")
286
287 }
288 else{
289
290 if (sw!=0 & swabs==0)
291 swabs <- as.integer(sw*nvoice)
292 if (swabs!=0 & sw==0)
293 sw <- swabs/nvoice
294
295 sllist <- checkarealsiglevel(sw,tw,w0,arealsiglevel,siglevel,1)
296 arealsiglevel <- sllist$arealsiglevel
297 siglevel <- sllist$siglevel
298
299 if (sw==0 & tw==0 & swabs==0) {
300 cat("# coherence without averaging makes no sense! \n")
301 siglevel <- 0
302 arealsiglevel <- 0
303 }
304
305 if (phase==FALSE) phs <- NULL
306
307 at <- NULL
308
309 tsadjust <- adjust.timeseries(ts1,ts2)
310 ts1 <- tsadjust$ts1
311 ts2 <- tsadjust$ts2
312
313 t <- time(ts1)
314 dt <- deltat(ts1)
315
316 s0rem <- s0
317 s0 <- adjust.s0(s0,dt)
318 dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1
319
320 noctave <- adjust.noctave(length(ts1),dt,s0,tw,noctave)
321
322 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
323
324 ts1anom <- ts1-mean(ts1)
325 ts2anom <- ts2-mean(ts2)
326
327 ts1.cwt <- cwt.ts(ts1anom,s0,noctave,nvoice,w0)
328 ts2.cwt <- cwt.ts(ts2anom,s0,noctave,nvoice,w0)
329
330 cosp <- Re(ts1.cwt)*Re(ts2.cwt) + Im(ts1.cwt)*Im(ts2.cwt)
331 quad <- Im(ts1.cwt)*Re(ts2.cwt) - Re(ts1.cwt)*Im(ts2.cwt)
332 wsp1 <- Mod(ts1.cwt)^2
333 wsp2 <- Mod(ts2.cwt)^2
334
335 smcosp <- smooth.matrix(cosp,swabs)
336 smquad <- smooth.matrix(quad,swabs)
337 smwsp1 <- smooth.matrix(wsp1,swabs)
338 smwsp2 <- smooth.matrix(wsp2,swabs)
339
340 smsmcosp <- smooth.time(smcosp,tw,dt,scalevector)
341 smsmquad <- smooth.time(smquad,tw,dt,scalevector)
342 smsmwsp1 <- smooth.time(smwsp1,tw,dt,scalevector)
343 smsmwsp2 <- smooth.time(smwsp2,tw,dt,scalevector)
344
345 if (sqr==FALSE)
346 wcoh <- sqrt((smsmcosp^2+smsmquad^2)/(smsmwsp1*smsmwsp2))
347 else
348 wcoh <- (smsmcosp^2+smsmquad^2)/(smsmwsp1*smsmwsp2)
349
350 if (phase)
351 phs <- atan2(smsmquad,smsmcosp)
352 else phs <- NULL
353
354 #POINTWISE SIGNIFICANCE TEST
355 if (siglevel[1]!=0) critval <- criticalvaluesWCO(s0,noctave,nvoice,w0,swabs,tw,siglevel)
356 else critval <- NULL
357 if (sqr==TRUE & is.null(critval)==FALSE)
358 critval <- critval^2
359
360 #AREAWISE SIGNIFICANCE TEST
361 if (arealsiglevel!=0){
362 atest <- arealtest(wcoh/critval,dt,s0,noctave,nvoice,w0,swabs,tw,siglevel,arealsiglevel,kernel,1)
363 at <- atest$at
364 kernel <- atest$kernel
365 }
366
367 wcoh[1,1] <- 0
368 wcoh[1,2] <- 1
369
370 if (phase){
371 phs[1,1] <- -pi
372 phs[1,2] <- pi
373 }
374
375 if (s0rem<s0){
376 wcoh <- addvalues(nvoice,dnoctave,wcoh,NA)
377 phs <- addvalues(nvoice,dnoctave,phs,NA)
378 noctave <- noctave+dnoctave
379 s0 <- s0/2^dnoctave
380 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
381 }
382
383 wclist <- list(modulus=wcoh,phase=phs,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=t,scales=scalevector,critval=critval,at=at,kernel=kernel)
384
385 class(wclist) <- "wt"
386
387 if (plot) plot(wclist,markt,marks,NULL,NULL,FALSE,phase,units,"Wavelet Coherence",device,file,split,color,pwidth,pheight,labsc,labtext,sigplot)
388
389 wclist
390
391 }
392
393 }
394
395#####################################
396# WCS
397#####################################
398
399wcs <- function(ts1,ts2,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,markt=-999,marks=-999,logscale=FALSE,phase=TRUE,plot=TRUE,units="",device="screen",file="wcsp",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext=""){
400
401 if (class(ts1)!="ts"){
402
403 cat("# This function needs two time series objects as input. You may construct them by using the function ts(data,start,deltat). Try '?ts' for help. \n")
404
405 }
406 else{
407
408 if (sw!=0 & swabs==0)
409 swabs <- as.integer(sw*nvoice)
410 if (swabs!=0 & sw==0)
411 sw <- swabs/nvoice
412
413 tsadjust <- adjust.timeseries(ts1,ts2)
414 ts1 <- tsadjust$ts1
415 ts2 <- tsadjust$ts2
416
417 t <- time(ts1)
418 dt <- deltat(ts1)
419
420 s0rem <- s0
421 s0 <- adjust.s0(s0,dt)
422 dnoctave <- as.integer(log(s0/s0rem-0.000000000001)/log(2))+1
423
424 noctave <- adjust.noctave(length(ts1),dt,s0,tw,noctave)
425
426 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
427
428 ts1anom <- ts1-mean(ts1)
429 ts2anom <- ts2-mean(ts2)
430
431 ts1.cwt <- cwt.ts(ts1anom,s0,noctave,nvoice,w0)
432 ts2.cwt <- cwt.ts(ts2anom,s0,noctave,nvoice,w0)
433
434 cosp <- Re(ts1.cwt)*Re(ts2.cwt) + Im(ts1.cwt)*Im(ts2.cwt)
435 quad <- Im(ts1.cwt)*Re(ts2.cwt) - Re(ts1.cwt)*Im(ts2.cwt)
436
437 smcosp <- smooth.matrix(cosp,swabs)
438 smquad <- smooth.matrix(quad,swabs)
439
440 smsmcosp <- smooth.time(smcosp,tw,dt,scalevector)
441 smsmquad <- smooth.time(smquad,tw,dt,scalevector)
442
443 wcsp <- smsmcosp^2+smsmquad^2
444
445 if (phase)
446 phs <- atan2(smsmquad,smsmcosp)
447 else phs <- NULL
448
449 if (phase){
450 phs[1,1] <- -pi
451 phs[1,2] <- pi
452 }
453
454 if (s0rem<s0){
455 wcsp <- addvalues(nvoice,dnoctave,wcoh,NA)
456 phs <- addvalues(nvoice,dnoctave,phs,NA)
457 noctave <- noctave+dnoctave
458 s0 <- s0/2^dnoctave
459 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
460 }
461
462 wclist <- list(modulus=wcsp,phase=phs,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=t,scales=scalevector,critval=NULL,at=NULL,kernel=NULL)
463
464 class(wclist) <- "wt"
465
466 if (plot) plot(wclist,markt,marks,NULL,NULL,logscale,phase,units,"Wavelet Cross Spectrum",device,file,split,color,pwidth,pheight,labsc,labtext)
467
468 wclist
469
470 }
471
472}
473
474##########################################
475# POINTWISE SIGNIFICANCE TEST
476##########################################
477
478rawWSP <- function(ts,s0=1,noctave=5,nvoice=20,w0=2*pi,swabs=0,tw=0,scalevector){
479
480 tsanom <- ts-mean(ts)
481
482 ts.cwt <- cwt.ts(tsanom,s0,noctave,nvoice,w0)
483
484 wsp <- Mod(ts.cwt)^2
485
486 smwsp <- smooth.matrix(wsp,swabs)
487 smsmwsp <- smooth.time(smwsp,tw,deltat(ts),scalevector)
488
489 smsmwsp
490
491}
492
493criticalvaluesWSP <- function(ts,s0=1,noctave=5,nvoice=10,w0=2*pi,swabs=0,tw=0,siglevel=0.95,nreal=1000){
494
495 t=time(ts)
496 dt=deltat(ts)
497
498 s0 <- adjust.s0(s0,dt)
499 noctave <- adjust.noctave(length(ts),dt,s0,tw,noctave)
500
501 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
502
503 cat("# calculating critical values by means of MC-simulations \n")
504
505 s0unit=s0/dt*w0/(2*pi)
506 s0log=as.integer((log2(s0unit)-1)*nvoice+1.5)
507
508 wsp0 <- rawWSP(ts,s0,noctave,nvoice,w0,swabs,tw,scalevector)
509
510 S <- dim(wsp0)[2]
511
512 n1 <- 1 + 2*as.integer( sqrt(2) * 2^((S-swabs+s0log-1)/nvoice+1) )
513 n2 <- max(scalevector)*tw*2/dt*1.1
514 n3 <- 2*tw*s0*2^noctave/dt+100
515 n <- max(n1,n2,n3)
516
517 center <- (n-1)/2+1
518 cv <- matrix(rep(0,S*length(siglevel)),ncol=S)
519 rmatrix <- matrix(0,nrow=nreal,ncol=S)
520
521 # Fitting AR1-process (red noise) to data
522 arts0 <- ar(ts,order.max=1)
523 sdts0 <- sd(ts[1:length(ts)])
524
525 if (arts0$order==0){
526 se <- sqrt(arts0$var)
527 arts0$ar <- 0.000001
528 }
529 else
530 se <- sqrt(sdts0*sdts0*(1-arts0$ar*arts0$ar))
531
532 tsMC <- ts(data=rep(0,n),start=t[1],frequency=1/dt)
533
534 # MC Simulations
535 for (i in 1:nreal){
536
537 tsMC[1:n] <- arima.sim(list(ar = arts0$ar), n+1, sd = se)[2:(n+1)]
538
539 rmatrix[i,] <- rawWSP(tsMC,s0,noctave,nvoice,w0,swabs,tw,scalevector)[center,]
540
541 }
542
543 for (s in (1+swabs):(S-swabs)) rmatrix[,s] <- sort(rmatrix[,s])
544
545 for (i in 1:length(siglevel)){
546 sigindex <- as.integer(nreal*siglevel[i])
547 cvv <- rmatrix[sigindex,]
548 cvv[is.na(cvv)] <- 0
549 cv[i,] <- cvv
550 }
551
552 cv
553
554}
555
556###########
557
558criticalvaluesWCO <- function(s0,noctave,nvoice,w0,swabs,tw,siglevel=0.95){
559
560 cv=rep(0,length(siglevel))
561
562 for (i in 1:length(siglevel)){
563
564 if (siglevel[i]!=0.9 && siglevel[i]!=0.95 && siglevel[i]!=0.99) siglevel[i] <- 0.95
565
566 if (siglevel[i]==0.9){
567 cat("# significance level set to 0.90 \n")
568 sw <- 1.0*swabs/nvoice
569 cv[i] <- 0.00246872*w0^2*sw + 0.0302419*w0*sw^2 + 0.342056*sw^3 -
570 0.000425853*w0^2 - 0.101749*w0*sw - 0.703537*sw^2 +
571 0.00816219*w0 + 0.442342*sw + 0.970315
572 }
573
574 if (siglevel[i]==0.95){
575 cat("# significance level set to 0.95 \n")
576 sw <- swabs*100.0/3.0/nvoice
577 cv[i] <- 0.0000823*w0^3 + 0.0000424*w0^2*sw + 0.0000113*w0*sw^2 +
578 0.0000154*sw^3 - 0.0023*w0^2 - 0.00219*w0*sw - 0.000751*sw^2 +
579 0.0205*w0 + 0.0127*sw + 0.95
580 }
581
582 if (siglevel[i]==0.99){
583 cat("# significance level set to 0.99 \n")
584 sw <- 1.0*swabs/nvoice
585 cv[i] <- 0.00102304*w0^2*sw + 0.00745772*w0*sw^2 + 0.230035*sw^3 -
586 0.000361565*w0^2 - 0.0502861*w0*sw - 0.440777*sw^2 +
587 0.00694998*w0 + 0.29643*sw + 0.972964
588 }
589
590 if (cv[i]>1) cv[i] <- 1
591
592 cat(paste("# significance testing, cv=",cv[i]," \n",sep=""))
593
594 }
595
596 cv
597
598}
599
600#############################
601# AREAWISE SIGNIFICANCE TEST
602#############################
603
604slide <- function(data,kernellist,s0,noctave,nvoice,cv){
605
606 # slides kernel over data set
607 #----------------------------
608 # data: matrix containing data
609 # kernellist: matrix containing kernel
610 # s0: lowest scale
611 # noctave: number of octaves
612 # nvoice: number of voices per octave
613 # cv: critical value, all values higher are set to one
614
615 #Time: (rows) n,i the kernel is to be scaled in this direction
616 #Scale: (cols) m,j
617
618 data[data<cv] <- 0
619
620 kernel <- kernellist$bitmap
621
622 js <- kernellist$is
623
624 sm <- s0*2^noctave
625
626 dbin <- tobin(data)
627 kbin <- tobin(kernel)
628
629 dn <- nrow(dbin)
630 dm <- ncol(dbin)
631
632 kn <- nrow(kbin)
633 km <- ncol(kbin)
634
635 mark <- matrix(rep(0,dn*dm),nrow=dn)
636
637 for (j in 1:(dm-km+1)){
638
639 s <- s0*2^((j+js-1)/nvoice)
640 kscn <- as.integer(kn*s/sm);
641 if (kscn==0) kscn <- 1
642
643 ksc <- scaleKernel(kbin,kscn)
644 kscm <- km
645
646 for (i in 1:(dn-kscn+1)){
647
648 subbin <- dbin[i:(i+kscn-1),j:(j+kscm-1)]
649
650 if (sum(ksc*subbin)==sum(ksc))
651 mark[i:(i+kscn-1),j:(j+kscm-1)] <- mark[i:(i+kscn-1),j:(j+kscm-1)]+ksc
652
653 }
654
655 }
656
657 mark <- tobin(mark)
658
659 mark
660
661}
662
663arealtest <- function(wt,dt=1,s0=1,noctave=5,nvoice=20,w0=2*pi,swabs=0,tw=0,siglevel,arealsiglevel=0.9,kernel=0,type=0){
664
665 slp <- slope(w0,swabs,tw,nvoice,siglevel,arealsiglevel,type)
666
667 if (length(kernel)<2){
668 maxarea <- s0*2^noctave*slp/10*nvoice/dt
669 cat(paste("# calculating kernel (maxarea=",maxarea,")\n",sep=""))
670 cvkernel <-
671 kernelRoot(s0,w0,a=maxarea,noctave,nvoice,swabs,tw,dt)
672
673 cat("# building kernel bitmap \n")
674 kernel <-
675 kernelBitmap(cvkernel,s0,w0,noctave,nvoice,swabs,tw,dt)
676
677 }
678
679 cat("# performing arealtest \n")
680 sl <- slide(wt,kernel,s0,noctave,nvoice,1)
681
682 list(at=sl,kernel=kernel)
683
684}
685
686#############################
687# PLOTTING
688#############################
689
690plotat <- function(t,wt,at,sigplot){
691
692 if (length(at)>1){
693 linewidth <- 1
694 if (sigplot==3)
695 linewidth <- 5
696
697 contour(x=t,z=at,levels=0.5,add=TRUE,drawlabels=FALSE,axes=FALSE,lwd=linewidth,col="black")
698 }
699
700}
701
702
703plotcv <- function(t,wt,cv){
704
705 if (length(dim(cv))==0)
706
707 contour(x=t,z=wt,levels=c(cv),drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1)
708
709 else{
710
711 for(i in 1:nrow(cv)){
712
713 v <- cv[i,]
714 v[v==0] <- 9999
715 m <- matrix(rep(v,length(t)),nrow=length(t),byrow=TRUE)
716
717 contour(x=t,z=wt/m,levels=1,drawlabels=FALSE,axes=FALSE,add=TRUE,col="black",lwd=1)
718
719 }
720
721 }
722
723}
724
725plotcoi <- function(t,s0,noctave,w0){
726
727 tv <- as.vector(t)
728 tvl <- tv[tv-tv[1]<(tv[length(tv)]-tv[1])/2]
729 tvr <- tv[tv-tv[1]>=(tv[length(tv)]-tv[1])/2]
730
731 lines(tvl,log2(((tvl-tv[1])*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black")
732 lines(tvr,log2(((tv[length(tv)]-tvr)*4*pi/((w0+sqrt(2+w0*w0))*sqrt(2)))/s0)/noctave,col="black")
733
734}
735
736plotmarks <- function(t,s0,noctave,markt,marks){
737
738 if (markt[1]!=-999){
739
740 for (i in 1:length(markt)){
741 lines(c(markt[i],markt[i]),c(0,1),lty="dotted")
742 }
743
744 }
745
746 if (marks[1]!=-999){
747
748 for (i in 1:length(marks)){
749 lines(c(t[1],t[length(t)]),c(log2(marks[i]/s0)/noctave,log2(marks[i]/s0)/noctave),lty="dotted")
750 }
751
752 }
753
754}
755
756
757#####################
758# PLOT.WT
759#####################
760
761plot.wt <- function(wt,markt=-999,marks=-999,t1=NULL,t2=NULL,logscale=FALSE,phase=FALSE,units="",plottitle="",device="screen",file="wt",split=FALSE,color=TRUE,pwidth=10,pheight=5,labsc=1.5,labtext="",sigplot=3,xax=NULL,xlab=NULL,yax=NULL,ylab=NULL){
762
763 plotwt(wt$modulus,wt$phase,wt$time,wt$s0,wt$noctave,wt$w0,wt$critval,wt$at,markt,marks,t1,t2,logscale,phase,units,plottitle,device,file,split,color,pwidth,pheight,labsc,labtext,sigplot,xax,xlab,yax,ylab)
764
765}
766
767plotwt <-
768 function(wt,phs,t,s0,noctave,w0,cv=NULL,at=NULL,markt=-999,marks=-999,t1=NULL,t2=NULL,logscale=FALSE,phase=FALSE,units="",plottitle="Wavelet Plot",device="screen",file="wavelet",split=FALSE,color=TRUE,pwidth=10,pheight=7,labsc=1,labtext="",sigplot=1,xax=NULL,xlab=NULL,yax=NULL,ylab=NULL){
769
770 if (is.null(phs)) phase <- FALSE
771
772 mgpv <- c(3,1,0)
773 if (labsc>1) mgpv[1] <- 3-(labsc-1.5)
774
775 ncolors <- 64
776 colors <- wtcolors(ncolors)
777 if (color==FALSE) colors <- gray((0:ncolors)/ncolors*0.7+0.3)
778
779 rangev <- (0:(ncolors-1))/(ncolors-1)
780 rangebar <- matrix(rangev,nrow=2,ncol=64,byrow=TRUE)
781
782 s <- 2^(0:(noctave))*s0
783 sn <- (0:(noctave))/noctave
784
785 deltat <- (t[length(t)]-t[1])/(length(t)-1)
786 cut <- FALSE
787 if ((!is.null(t1)) | (!is.null(t2))){
788 if (t1<t2 & t1>=t[1] & t2<=t[length(t)]){
789
790 cut <- TRUE
791
792 i1 <- (t1-t[1])/deltat+1
793 i2 <- (t2-t[1])/deltat+1
794
795 t <- t[t>=t1 & t<=t2]
796
797 wt <- wt[i1:i2,]
798 if (phase) phs <- phs[i1:i2,]
799 if (length(at)>1) at <- at[i1:i2,]
800
801 }
802 }
803
804 #----------------
805 # Setting layout
806 #----------------
807
808 if (device=="ps" && split==FALSE){
809
810 file <- paste(file,".eps",sep="")
811
812 postscript(file,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight)
813 cat(paste("# writing",file," \n"))
814
815 }
816
817 if (device=="ps" && split==TRUE){
818
819 file1 <- paste(file,".mod.eps",sep="")
820
821 postscript(file1,onefile=TRUE,horizontal=TRUE,paper="special",width=pwidth,height=pheight)
822 cat(paste("# writing",file1," \n"))
823
824 }
825
826 if (phase==TRUE && split==FALSE) nlo <- layout(matrix(c(1,2,3,4),2,2,byrow=TRUE),widths=c(4,1))
827 else nlo <- layout(matrix(c(1,2),ncol=2,byrow=TRUE),widths=c(4,1))
828
829
830 #-----------------
831 # Plotting modulus
832 #-----------------
833
834 if (logscale){
835 if (units=="")
836 image(x=t,z=log10(wt),col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
837 else
838 image(x=t,z=log10(wt),col = colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
839 }
840 else{
841 if (units=="")
842 image(x=t,z=wt,col = colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
843 else
844 image(x=t,z=wt,col = colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
845 }
846
847 text(t[1+as.integer(length(t)*0.1)],0.8,labtext,cex=2)
848
849 if (sigplot==1 | sigplot==3){
850 if (is.null(cv)==FALSE){ #Critical values
851 if (cv[1]!=0 & cv[1]!=-1) plotcv(t,wt,cv)
852 }
853 }
854
855 if (sigplot>1) plotat(t,wt,at,sigplot)
856
857 if (!cut) plotcoi(t,s0,noctave,w0) #cone of influence
858 plotmarks(t,s0,noctave,markt,marks) #additional guiding lines
859
860 if (is.null(xax))
861 axis(side=1,at=axTicks(1),cex.axis=labsc)
862 else
863 if (is.null(xlab))
864 axis(side=1,xax,labels=as.character(xax),cex.axis=labsc)
865 else
866 axis(side=1,xax,labels=xlab,cex.axis=labsc)
867
868 if (is.null(yax))
869 axis(side=2,sn,labels=as.character(s),cex.axis=labsc)
870 else
871 if (is.null(ylab))
872 axis(side=2,yax,labels=as.character(yax),cex.axis=labsc)
873 else
874 axis(side=2,yax,labels=ylab,cex.axis=labsc)
875
876 title(main=plottitle,cex=labsc)
877
878 image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
879
880 if (is.null(cv)==FALSE){
881 if (length(dim(cv))==0){
882 for (i in 1:length(cv))
883 if (cv[i]!=0) lines(c(-1,2),c(cv[i],cv[i]))
884 }
885 }
886
887 if (!logscale)
888 axis(side=2,(0:5)/5,labels=c("0","","","","","1"),cex.axis=labsc)
889 else{
890 labelv <- substr(as.character(c(0:5)*(max(log10(wt),na.rm=TRUE)-min(log10(wt),na.rm=TRUE))/5),1,4)
891 axis(side=2,(0:5)/5,labels=labelv,cex.axis=labsc)
892 }
893
894
895 #-----------------
896 # Plotting phase
897 #-----------------
898 if (phase==TRUE){
899
900 if (device=="ps" && split==TRUE){
901
902 dev.off()
903
904 file2 <- paste(file,".phs.eps",sep="")
905
906 postscript(file2,onefile=TRUE,horizontal=TRUE,paper="special",width=10,height=5)
907 cat(paste("# writing",file2," \n"))
908
909 }
910
911 colors <- rainbow(ncolors)
912 if (color==FALSE){
913 dummy <- gray((0:ncolors)/ncolors)
914 colors[1:(ncolors/2)] <- dummy[(ncolors/2+1):ncolors]
915 colors[(ncolors/2+1):ncolors] <- dummy[1:(ncolors/2)]
916 }
917
918 if (units=="")
919 image(x=t,z=phs,col=colors,axes=FALSE,xlab="Time",ylab="Scale",frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
920 else
921 image(x=t,z=phs,col=colors,axes=FALSE,xlab=paste("Time ","[",units,"]",sep=""),ylab=paste("Scale ","[",units,"]",sep=""),frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
922
923 if (is.null(cv)==FALSE) plotcv(t,wt,cv)
924 if (sigplot>1) plotat(t,wt,at,sigplot)
925 if (!cut) plotcoi(t,s0,noctave,w0)
926 plotmarks(t,s0,noctave,markt,marks)
927
928 if (is.null(xax))
929 axis(side=1,at=axTicks(1),cex.axis=labsc)
930 else
931 if (is.null(xlab))
932 axis(side=1,xax,labels=as.character(xax),cex.axis=labsc)
933 else
934 axis(side=1,xax,labels=xlab,cex.axis=labsc)
935
936 if (is.null(yax))
937 axis(side=2,sn,labels=as.character(s),cex.axis=labsc)
938 else
939 if (is.null(ylab))
940 axis(side=2,yax,labels=as.character(yax),cex.axis=labsc)
941 else
942 axis(side=2,yax,labels=ylab,cex.axis=labsc)
943
944
945 title(main="Phase")
946
947 image(z=rangebar,axes=FALSE,col=colors,frame.plot=TRUE,cex.lab=labsc,mgp=mgpv)
948 axis(side=2,(0:4)/4,labels=c("-PI","","0","","PI"),cex.axis=labsc)
949
950 }
951
952 if (device=="ps") dev.off()
953
954 }
955
956##############################
957# Surrogates
958##############################
959
960createwavesurrogates <- function(nsurr=1,wt=1,n,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){
961
962 surrmat <- matrix(rep(0,n*nsurr),ncol=nsurr)
963
964 for (i in 1:nsurr){
965
966 cat(paste("# Creating realization ",i,"\n",sep=""))
967
968 x <- rnorm(n)
969 xts <- ts(x,start=0,deltat=dt)
970
971 xwt <- cwt.ts(xts,s0,noctave,nvoice,w0)
972 wtsurr <- wt*xwt
973
974 surri <- as.vector(invmorlet(wtsurr,0,dt,s0,noctave,nvoice,w0))
975
976 surrmat[,i] <- Re(surri)
977
978 }
979
980 surrmat
981
982}
983
984surrwave <- function(x,...)
985 UseMethod("surrwave")
986
987surrwave.wt <- function(wt,nsurr=1,spec=TRUE){
988
989 n <- length(wt$time)
990 t0 <- wt$time[1]
991 dt <- (wt$time[13]-t0)/12
992 s0 <- wt$s0
993 noctave <- wt$noctave
994 nvoice <- wt$nvoice
995 w0 <- wt$w0
996
997 wt <- wt$modulus
998 if (spec==TRUE)
999 wt <- sqrt(wt)
1000
1001 surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0)
1002
1003 surrts <- ts(surrmat,start=t0,deltat=dt)
1004
1005 surrts
1006
1007}
1008
1009surrwave.matrix <- function(mat,nsurr=1,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,spec=FALSE){
1010
1011 if (sw!=0 & swabs==0)
1012 swabs <- as.integer(sw*nvoice)
1013
1014 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
1015
1016 if ((noctave*nvoice+1)!=dim(mat)[2])
1017 cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n")
1018
1019 n <- dim(mat)[1]
1020
1021 if (spec==FALSE)
1022 mat <- Mod(mat)
1023 else
1024 mat <- sqrt(Mod(mat))
1025
1026 wtsm <- smooth.matrix(mat,swabs)
1027 wt <- smooth.time(wtsm,tw,dt,scalevector)
1028
1029 surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0)
1030
1031 surrts <- ts(surrmat,start=t0,deltat=dt)
1032
1033 surrts
1034
1035}
1036
1037surrwave.character <-
1038 function(file,nsurr=1,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0,transpose=TRUE,spec=FALSE){
1039
1040 if (sw!=0 & swabs==0)
1041 swabs <- as.integer(sw*nvoice)
1042
1043 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
1044
1045 if (transpose==FALSE)
1046 mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=TRUE)
1047 else
1048 mat <- matrix(scan(file,comment.char="#"),ncol=nvoice*noctave+1,byrow=FALSE)
1049
1050 if ((noctave*nvoice+1)!=dim(mat)[2])
1051 cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n")
1052
1053 n <- dim(mat)[1]
1054
1055 if (spec==FALSE)
1056 mat <- Mod(mat)
1057 else
1058 mat <- sqrt(Mod(mat))
1059
1060 wtsm <- smooth.matrix(mat,swabs)
1061 wt <- smooth.time(wtsm,tw,dt,scalevector)
1062
1063 surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0)
1064
1065 surrts <- ts(surrmat,start=t0,deltat=dt)
1066
1067 surrts
1068
1069 }
1070
1071surrwave.ts <- function(ts,nsurr=1,s0=1,noctave=5,nvoice=10,w0=2*pi,sw=0,tw=0,swabs=0){
1072
1073 n <- length(ts)
1074 t0 <- time(ts)[1]
1075 dt <- deltat(ts)
1076 if (sw!=0 & swabs==0)
1077 swabs <- as.integer(sw*nvoice)
1078
1079 scalevector <- 2^(0:(noctave*nvoice)/nvoice)*s0
1080
1081 wt <- Mod(cwt.ts(ts,s0,noctave,nvoice,w0))
1082
1083 wtsm <- smooth.matrix(wt,swabs)
1084 wt <- smooth.time(wtsm,tw,dt,scalevector)
1085
1086 surrmat <- createwavesurrogates(nsurr,wt,n,dt,s0,noctave,nvoice,w0)
1087
1088 surrts <- ts(surrmat,start=t0,deltat=dt)
1089
1090 surrts
1091
1092}
1093
1094invmorlet <- function(wt,t0=0,dt=1,s0=1,noctave=5,nvoice=10,w0=2*pi){
1095
1096 if ((noctave*nvoice+1)!=dim(wt)[2])
1097 cat("# ERROR! nscales unequal noctave*nvoice+1 ! \n")
1098
1099 n <- dim(wt)[1]
1100
1101 wt[is.na(wt)] <- 0
1102
1103 tsRe <- rep(0,n)
1104 tsIm <- rep(0,n)
1105
1106 wtRe <- t(Re(wt))
1107 wtIm <- t(Im(wt))
1108
1109 z <- .C("invmorlet",
1110 as.double(as.vector(wtRe)),
1111 as.double(as.vector(wtIm)),
1112 as.integer(n),
1113 as.double(dt),
1114 as.double(s0),
1115 as.integer(noctave),
1116 as.integer(nvoice),
1117 as.double(w0),
1118 tsRetmp = as.double(tsRe),
1119 tsImtmp = as.double(tsIm),
1120 PACKAGE="sowas")
1121
1122 invvec=complex(real=z$tsRetmp,imaginary=z$tsImtmp)
1123 invts <- ts(data=invvec,start=t0,deltat=dt)
1124
1125 invts
1126
1127}
1128
1129#################################
1130# INPUT / OUTPUT
1131#################################
1132
1133readmatrix <- function(file,M){
1134
1135 A <- matrix(scan(file,comment.char="#"),ncol=M,byrow=TRUE)
1136
1137 A
1138
1139}
1140
1141
1142readts <- function(file){
1143
1144 A <- matrix(scan(file,comment.char="#"),ncol=2,byrow=TRUE)
1145
1146 Adum <- A
1147
1148 Adum[is.na(Adum)] <- 0
1149
1150 t <- Adum %*% c(1,0)
1151 x <- A %*% c(0,1)
1152
1153 N=length(t)
1154 f=1/(t[13]-t[1])*12
1155
1156 if ((f>11) && (f<13)) f <- 12
1157
1158 timeseries<-ts(data=x,start=t[1],frequency=f)
1159
1160 timeseries
1161
1162}
1163
1164writematrix <- function(file,data,header="# R Matrix"){
1165
1166 write(header,file)
1167 write(t(data),file,ncol(data),append=TRUE)
1168
1169}
1170
1171############################
1172# HELP FUNCTIONS
1173############################
1174
1175smooth.matrix <- function(wt,swabs){
1176
1177 if (swabs != 0)
1178 smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1)))
1179 else
1180 smwt <- wt
1181
1182 smwt
1183
1184}
1185
1186smooth.time <- function(wt,tw,dt,scalevector){
1187
1188 smwt <- wt
1189
1190 if (tw != 0){
1191 for (i in 1:length(scalevector)){
1192
1193 twi <- as.integer(scalevector[i]*tw/dt)
1194 smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1))
1195
1196 }
1197 }
1198
1199 smwt
1200
1201}
1202
1203
1204adjust.noctave <- function(N,dt,s0,tw,noctave){
1205
1206 if (tw>0){
1207 dumno <- as.integer((log(N*dt)-log(2*tw*s0))/log(2))
1208 if (dumno<noctave){
1209 cat("# noctave adjusted to time smoothing window \n")
1210 noctave <- dumno
1211 }
1212 }
1213
1214 noctave
1215
1216}
1217
1218adjust.s0 <- function(s0,dt){
1219
1220 if (s0<2*dt){
1221 s0 <- 2*dt
1222 cat(paste("# s0 set to ",s0," \n"))
1223 }
1224
1225 s0
1226
1227}
1228
1229adjust.timeseries <- function(ts1,ts2){
1230
1231 if (length(ts1)!=length(ts2)){
1232 tsint <- ts.intersect(ts1,ts2)
1233 dt <- deltat(ts1)
1234 ts1 <- ts(data=tsint[,1],start=time(tsint)[1],frequency=1/dt)
1235 ts2 <- ts(data=tsint[,2],start=time(tsint)[1],frequency=1/dt)
1236 t <- time(ts1)
1237 }
1238
1239 list(ts1=ts1,ts2=ts2)
1240
1241}
1242
1243checkarealsiglevel <- function(sw,tw,w0,arealsiglevel,siglevel,type){
1244
1245 if (type==0){
1246
1247 swv <- c(0,0.5,1)
1248 twv <- c(0,1.5,3)
1249 w0v <- c(pi,2*pi)
1250
1251 if (length(swv[swv==sw])==0 || length(twv[twv==tw])==0 ||
1252 length(w0v[w0v==w0])==0){
1253 arealsiglevel <- 0
1254 cat("# areawise test for spectrum currently \n")
1255 cat("# only possible for \n")
1256 cat("# sw = 0 \n")
1257 cat("# tw = 0 \n")
1258 cat("# w0 = 2pi \n")
1259 cat("# No areawise test performed \n")
1260 }
1261 }
1262
1263 if (type==1){
1264
1265 swv <- c(0.5)
1266 twv <- c(1.5)
1267 w0v <- c(2*pi)
1268
1269 if (length(swv[swv==sw])==0 || length(twv[twv==tw])==0 ||
1270 length(w0v[w0v==w0])==0){
1271 arealsiglevel <- 0
1272 cat("# areawise test for coherence currently \n")
1273 cat("# only possible for \n")
1274 cat("# sw = 0.5 \n")
1275 cat("# tw = 1.5 \n")
1276 cat("# w0 = 2pi \n")
1277 cat("# No areawise test performed \n")
1278 }
1279 }
1280
1281
1282 if (arealsiglevel!=0){
1283 arealsiglevel <- 0.9
1284 siglevel <- 0.95
1285 cat("# currently only siglevel=0.95 and arealsiglevel=0.9 possible for areawise test \n")
1286 }
1287
1288 list(siglevel=siglevel,arealsiglevel=arealsiglevel)
1289
1290}
1291
1292########################
1293
1294as.wt <- function(modulus,phase=NULL,s0=NULL,noctave=NULL,nvoice=NULL,w0=NULL,dt=1,time=NULL,scales=NULL,critval=NULL,at=NULL,kernel=NULL,N=NULL,t0=NULL){
1295
1296 if (is.null(scales))
1297 gotscales <- FALSE
1298 else
1299 gotscales <- TRUE
1300
1301 if ((!gotscales) & (!is.null(s0)) & (!is.null(noctave)) &(!is.null(nvoice))){
1302 gotscales <- TRUE
1303 scales=2^(0:(noctave*nvoice)/nvoice)*s0
1304 }
1305
1306 if (gotscales & (is.null(s0) | is.null(noctave) |
1307 is.null(nvoice))){
1308 s0 <- scales[1]
1309 noctave <- log(scales[length(scales)]/s0)/log(2)
1310 nvoice <- (length(scales)-1)/noctave
1311 }
1312
1313
1314 if (!gotscales)
1315 cat("# ERROR! No scales given! \n")
1316
1317 if (is.null(time)) gottimes <- FALSE
1318 else gottimes <- TRUE
1319
1320 if ((!gottimes) & (!is.null(dt)) & (!is.null(t0)) &(!is.null(N))){
1321 gottimes <- TRUE
1322 time=(0:(N-1))*dt+t0
1323 }
1324
1325 if (!gottimes)
1326 cat("# ERROR! No time vector given! \n")
1327
1328 wcolist <- list(modulus=modulus,phase=phase,s0=s0,noctave=noctave,nvoice=nvoice,w0=w0,time=time,scales=scales,critval=critval,at=at,kernel=kernel)
1329
1330 class(wcolist) <- "wt"
1331
1332 wcolist
1333
1334}
1335
1336########################
1337
1338wtcolors <- function(ncolors){
1339
1340 upside <- rainbow(ncolors,start=0,end=.7)
1341 #upside <- heat.colors(ncolors+5)
1342 #upside <- upside[1:ncolors]
1343
1344
1345 down <- upside
1346
1347 for (i in 1:ncolors){
1348 down[i] <- upside[ncolors+1-i]
1349 }#
1350
1351 down
1352
1353}
1354
1355####################
1356
1357createwgn <- function(N,sig,dt){
1358
1359 timeseries<-ts(data=rnorm(N,0,sig),start=0,deltat=dt)
1360
1361 timeseries
1362
1363}
1364
1365
1366createar <- function(N,a,sig,dt){
1367
1368 if (a==0) a <- 0.000000001
1369
1370 se <- sqrt(sig*sig*(1-a*a))
1371 tsMC <- ts(data=rep(0,N),start=0,deltat=dt)
1372
1373 tsMC[1:N] <- arima.sim(list(ar = a), N, sd = se)
1374
1375}
1376
1377######################
1378
1379rk <- function(N=1000,s=8,noctave=5,nvoice=10,w0=2*pi,plot=TRUE){
1380
1381 t <- 1:N
1382
1383 sunit <- s*(w0+sqrt(2+w0*w0))/(4*pi)
1384
1385 s0 <- 4
1386 #s0unit <- s0*(w0+sqrt(2+w0*w0))/(4*pi)
1387 s0unit=s0/dt*w0/(2*pi) #(CORRECT)
1388 s0log <- as.integer((log2(s0unit)-1)*nvoice+1.5)
1389
1390 totnoct <- noctave+as.integer(s0log/nvoice)+1
1391
1392 x <- morlet(N,N/2,sunit,w0)
1393
1394 totts.cwt <- Mod(cwt(x,totnoct,nvoice,w0,plot=0))
1395 wt=totts.cwt[,s0log:(s0log+noctave*nvoice)]
1396
1397 wt <- wt/max(wt)
1398
1399 if (plot==TRUE) plotwt(wt,0,t,s0,noctave,w0,units="",plottitle="Reproducing Kernel")
1400
1401 wt
1402
1403}
1404
1405###################
1406
1407addvalues <- function(nvoice,dnoctave,x,value){
1408
1409 nr <- dim(x)[1]
1410 nc <- dim(x)[2]
1411 dnc <- nvoice*dnoctave
1412
1413 y <- matrix(rep(value,nr*(nc+dnc)),nrow=nr)
1414
1415 y[,(dnc+1):(dnc+nc)] <- x
1416
1417 y
1418
1419}
1420
1421####################
1422
1423scalematrix <- function(wt){
1424
1425 # scales matrix, such that the maximal value is one
1426 # wt: matrix to be scaled
1427
1428 mval <- max(wt,na.rm=TRUE)
1429
1430 wt <- wt/mval
1431
1432 wt
1433
1434}
1435
1436
1437foldKernel <- function(F,swabs,tw,s,dt){
1438
1439 # folds a matrix (e.g. kernel with smoothing window
1440 # F: matrix input
1441 # swabs: smooting window width
1442
1443 smsF <- smooth.matrix(F,swabs)
1444 smsF[is.na(smsF)] <- 0
1445
1446 smtsF <- smooth.time(smsF,tw,dt,s)
1447
1448 smtsF[is.na(smtsF)] <- 0
1449
1450 scF <- scalematrix(smtsF)
1451
1452 scF
1453
1454}
1455
1456kernelBitmap <- function(c,s0=1,w0=6,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){
1457 # produces kernel bitmap
1458 # c: height of contour, that defines area
1459 # s0: lowest scale
1460 # noctave: number of octaves
1461 # nvoice: number of voices per octave
1462 # swabs: smoothing window length in scale direction
1463 # dt: sampling time
1464
1465 s <- s0*2^noctave
1466 is <- noctave*nvoice
1467
1468 x <- s0*2^(((1:(nvoice*(noctave+2)))-1)/nvoice)
1469
1470 t <- max(2000,max(x)*tw*2/dt*1.1)
1471
1472 y <- ((0:(2*t))-t)/2*dt
1473
1474 X <- matrix(x,ncol=nvoice*(noctave+2),nrow=2*t+1,byrow=T)
1475 Y <- matrix(y,ncol=nvoice*(noctave+2),nrow=2*t+1)
1476
1477 F <- sqrt(2*s*X/(s*s+X*X))*exp(-0.5*(Y*Y+w0*w0*(X-s)*(X-s))/(s*s+X*X))
1478
1479 F <- foldKernel(F,swabs,tw,x,dt)
1480
1481 F[F<c] <- 0
1482 F[F>=c] <- 1
1483
1484 is1 <- 1
1485 is2 <- nvoice*(noctave+1)
1486 it1 <- 1
1487 it2 <- 2*t+1
1488
1489 L <- F[1:(2*t+1),is1]
1490
1491 while (length(L[L!=0])==0) {
1492 is1 <- is1+1
1493 L <- F[1:(2*t+1),is1]
1494 }
1495
1496 L <- F[1:(2*t+1),is2]
1497
1498 while (length(L[L!=0])==0) {
1499 is2 <- is2-1
1500 L <- F[1:(2*t+1),is2]
1501 }
1502
1503
1504 L <- F[it1,1:(nvoice*(noctave+2))]
1505
1506 while (length(L[L!=0])==0) {
1507 it1 <- it1+1
1508 L <- F[it1,1:(nvoice*(noctave+2))]
1509 }
1510
1511 L <- F[it2,1:(nvoice*(noctave+2))]
1512
1513 while (length(L[L!=0])==0) {
1514 it2 <- it2-1
1515 L <- F[it2,1:(nvoice*(noctave+2))]
1516 }
1517
1518 kernel <- list(bitmap=F[(it1-1):(it2+1),(is1-1):(is2+1)],is=is-is1)
1519
1520 kernel
1521
1522}
1523
1524kernelRoot <- function(s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){
1525
1526 tol <- 0.005
1527 cmin <- 0
1528 cmax <- 1
1529 cntr <- 0.5
1530 da <- a
1531
1532 while (abs(da/a)>tol){
1533
1534 da <- kernelArea(cntr,s0,w0,a,noctave,nvoice,swabs,tw,dt)
1535 if (da>0){
1536 cmin <- cntr
1537 cntr <- mean(c(cntr,cmax))
1538 }
1539 if (da<0){
1540 cmax <- cntr
1541 cntr <- mean(c(cntr,cmin))
1542 }
1543 }
1544
1545 cntr
1546
1547}
1548
1549kernelArea <- function(cntr,s0=1,w0=6,a=0,noctave=6,nvoice=20,swabs=0,tw=0,dt=1){
1550
1551 # calulates area of reproducing kernel for smoothed data at scale s0*2^noctave
1552 # cntr: height of contour line to define kernel area. This
1553 # parameter is to be estimated!
1554 # s0: lowest scale
1555 # w0: parameter of Morlet Wavelet
1556 # a: area offset (only needed, when finding root. Is set to
1557 # desired area
1558 # noctave: number of octaves
1559 # nvoice: number of voices per octave
1560 # swabs: smoothing window width in scale direction
1561 # dt: sampling time
1562
1563 s <- s0*2^noctave
1564
1565 x <- s0*2^(((1:(nvoice*(noctave+2)))-1)/nvoice)
1566
1567 t <- max(2000,max(x)*tw*2/dt*1.1)
1568
1569 y <- ((0:(2*t))-t)/2*dt
1570
1571 X <- matrix(x,ncol=nvoice*(noctave+2),nrow=2*t+1,byrow=T)
1572 Y <- matrix(y,ncol=nvoice*(noctave+2),nrow=2*t+1)
1573
1574 F <- sqrt(2*s*X/(s*s+X*X))*exp(-0.5*(Y*Y+w0*w0*(X-s)*(X-s))/(s*s+X*X))
1575
1576 F <- foldKernel(F,swabs,tw,x,dt)
1577
1578 F[F>=cntr] <- 1
1579 F[F<cntr] <- 0
1580
1581 area <- length(F[F==1])-a
1582
1583 area
1584
1585}
1586
1587
1588tobin <- function(x){
1589
1590 # sets nonzero values to one
1591
1592 y <- x/x
1593 y[is.na(y)] <- 0
1594
1595 y
1596
1597}
1598
1599scaleKernel <- function(kernel,l){
1600
1601 # scales kernel length in time direction proportional to scale
1602 # kernel: data bitmap of width n
1603 # l: new width of kernel
1604
1605 n <- nrow(kernel)
1606 m <- ncol(kernel)
1607
1608 newKernel <- matrix(rep(0,m*l),nrow=l)
1609
1610 d <- as.double(n)/as.double(l)
1611
1612 for (i in 1:l){
1613 j <- as.integer((i-0.5)*d)
1614 if (j==0) j <- 1
1615 newKernel[i,1:m] <- kernel[j,1:m]
1616 }
1617
1618 newKernel
1619
1620}
1621
1622slope <- function(w0,swabs,tw,nvoice,siglevel,arealsiglevel,type){
1623
1624 sw <- swabs/nvoice
1625
1626 if (type==0){ # wavelet spectrum
1627
1628 if (tw == 0 & sw == 0 & w0 == 1 *pi) slp <- 5.82518 # w = 18.35831
1629 if (tw == 1.5 & sw == 0 & w0 == 1 *pi) slp <- 24.69852 # w = 14.30709
1630 if (tw == 3 & sw == 0 & w0 == 1 *pi) slp <- 35.48368 # w = 14.72354
1631 if (tw == 0 & sw == 5 & w0 == 1 *pi) slp <- 7.347707 # w = 17.96942
1632 if (tw == 1.5 & sw == 5 & w0 == 1 *pi) slp <- 28.24291 # w = 12.65993
1633 if (tw == 3 & sw == 5 & w0 == 1 *pi) slp <- 51.13723 # w = 10.96359
1634 if (tw == 0 & sw == 10 & w0 == 1 *pi) slp <- 10.47856 # w = 15.5941
1635 if (tw == 1.5 & sw == 10 & w0 == 1 *pi) slp <- 45.07387 # w = 15.29793
1636 if (tw == 3 & sw == 10 & w0 == 1 *pi) slp <- 52.82886 # w = 12.72361
1637
1638 if (tw == 0 & sw == 0 & w0 == 2 *pi) slp <- 8.718912 # w = 17.75933
1639 if (tw == 1.5 & sw == 0 & w0 == 2 *pi) slp <- 11.88006 # w = 15.39648
1640 if (tw == 3 & sw == 0 & w0 == 2 *pi) slp <- 26.55977 # w = 1.064384
1641 if (tw == 0 & sw == 5 & w0 == 2 *pi) slp <- 14.64761 # w = 16.27518
1642 if (tw == 1.5 & sw == 5 & w0 == 2 *pi) slp <- 28.27798 # w = 14.57059
1643 if (tw == 3 & sw == 5 & w0 == 2 *pi) slp <- 63.54121 # w = 12.83778
1644 if (tw == 0 & sw == 10 & w0 == 2 *pi) slp <- 27.78735 # w = 11.95813
1645 if (tw == 1.5 & sw == 10 & w0 == 2 *pi) slp <- 41.27260 # w = 12.03379
1646 if (tw == 3 & sw == 10 & w0 == 2 *pi) slp <- 67.37015 # w = 10.63935
1647
1648 }
1649
1650 if (type==1){ # wavelet coherence
1651
1652 if (tw==0 & sw==0 & w0==pi) slp <- 999 #not valid
1653 if (tw==1.5 & sw==0 & w0==pi) slp <- 1
1654 if (tw==3 & sw==0 & w0==pi) slp <- 1
1655 if (tw==0 & sw==0.5 & w0==pi) slp <- 1
1656 if (tw==1.5 & sw==0.5 & w0==pi) slp <- 1
1657 if (tw==3 & sw==0.5 & w0==pi) slp <- 1
1658 if (tw==0 & sw==1 & w0==pi) slp <- 1
1659 if (tw==1.5 & sw==1 & w0==pi) slp <- 1
1660 if (tw==3 & sw==1 & w0==pi) slp <- 1
1661
1662 if (tw==0 & sw==0 & w0==2*pi) slp <- 999 #not valid
1663 if (tw==1.5 & sw==0 & w0==2*pi) slp <- 1
1664 if (tw==3 & sw==0 & w0==2*pi) slp <- 1
1665 if (tw==0 & sw==0.5 & w0==2*pi) slp <- 1
1666 if (tw==1.5 & sw==0.5 & w0==2*pi) slp <- 8.3
1667 if (tw==3 & sw==0.5 & w0==2*pi) slp <- 1
1668 if (tw==0 & sw==1 & w0==2*pi) slp <- 1
1669 if (tw==1.5 & sw==1 & w0==2*pi) slp <- 1
1670 if (tw==3 & sw==1 & w0==2*pi) slp <- 1
1671
1672 if (tw==0 & sw==0 & w0==3*pi) slp <- 999 #not valid
1673 if (tw==1.5 & sw==0 & w0==3*pi) slp <- 1
1674 if (tw==3 & sw==0 & w0==3*pi) slp <- 1
1675 if (tw==0 & sw==0.5 & w0==3*pi) slp <- 1
1676 if (tw==1.5 & sw==0.5 & w0==3*pi) slp <- 1
1677 if (tw==3 & sw==0.5 & w0==3*pi) slp <- 1
1678 if (tw==0 & sw==1 & w0==3*pi) slp <- 1
1679 if (tw==1.5 & sw==1 & w0==3*pi) slp <- 1
1680 if (tw==3 & sw==1 & w0==3*pi) slp <- 1
1681
1682 if (tw==0 & sw==0 & w0==4*pi) slp <- 999 #not valid
1683 if (tw==1.5 & sw==0 & w0==4*pi) slp <- 1
1684 if (tw==3 & sw==0 & w0==4*pi) slp <- 1
1685 if (tw==0 & sw==0.5 & w0==4*pi) slp <- 1
1686 if (tw==1.5 & sw==0.5 & w0==4*pi) slp <- 1
1687 if (tw==3 & sw==0.5 & w0==4*pi) slp <- 1
1688 if (tw==0 & sw==1 & w0==4*pi) slp <- 1
1689 if (tw==1.5 & sw==1 & w0==4*pi) slp <- 1
1690 if (tw==3 & sw==1 & w0==4*pi) slp <- 1
1691
1692 }
1693
1694 cat(paste("# slope ",slp,"\n",sep=""))
1695
1696 slp
1697
1698}
1699
1700#============================
1701#============================
1702#============================
1703
1704## 1. Read auxiliar data files ####
1705
1706identifiants <- read.table("identifs.txt")[ ,1]
1707dates0 <- read.table("datesall.txt")[, 1]
1708dates <- as.character(dates0[grep("2009", dates0)])
1709rm(dates0)
1710
1711n <- length(identifiants)
1712p <- delta <- length(dates)
1713
1714synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200RC.txt")))
1715#synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200-random.txt")))
1716
1717nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing
1718synchros09[nas, 1] <- rowMeans(synchros09[nas, 2:4]) #valeurs après 1er janvier
1719
1720#moyenne pondérée pour compléter deux demi-heures manquantes
1721imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
1722synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518])
1723
1724conso <- synchros09[-201, ]; # series must be on rows
1725n <- nrow(conso)
1726delta <- ncol(conso)
1727
1728rm(synchros09, nas)
1729
1730## 2. Compute WER distance matrix ####
1731
1732## _.a CWT -- Filtering the lowest freqs (>6m) ####
1733# nvoice <- 4
1734# # noctave4 = 2^13 = 8192 half hours ~ 180 days
1735# noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2,
1736# tw = 0, noctave = 13)
1737# # 4 here represent 2^5 = 32 half-hours ~ 1 day
1738# scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2
1739# lscvect4 <- length(scalevector4)
1740# lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect
1741
1742
1743#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]
1744
1745#TODO: une fonction qui fait lignes 59 à 91
1746
1747#cube:
1748# Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
1749# scalevector = scalevector4,
1750# lt = delta, smooth = FALSE,
1751# nvoice = nvoice) # observations node with CWT
1752#
1753# #matrix:
1754# ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
1755# #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
1756#
1757# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
1758# for(i in 1:n)
1759# Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
1760#
1761# #rm(conso, Xcwt4); gc()
1762#
1763# ## _.b WER^2 distances ########
1764# Xwer_dist <- matrix(0.0, n, n)
1765# for(i in 1:(n - 1)){
1766# mat1 <- vect2mat(Xcwt2[i,])
1767# for(j in (i + 1):n){
1768# mat2 <- vect2mat(Xcwt2[j,])
1769# num <- Mod(mat1 * Conj(mat2))
1770# WX <- Mod(mat1 * Conj(mat1))
1771# WY <- Mod(mat2 * Conj(mat2))
1772# smsmnum <- smCWT(num, scalevector = scalevector4)
1773# smsmWX <- smCWT(WX, scalevector = scalevector4)
1774# smsmWY <- smCWT(WY, scalevector = scalevector4)
1775# wer2 <- sum(colSums(smsmnum)^2) /
1776# sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
1777# Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
1778# Xwer_dist[j, i] <- Xwer_dist[i, j]
1779# }
1780# }
1781# diag(Xwer_dist) <- numeric(n)
1782#
1783# save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
1784# save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")
1785
1786load("../res/2009_synchros200WER.Rdata")
1787#load("../res/2009_synchros200-randomWER.Rdata")
1788
1789## 3. Cluster using WER distance matrix ####
1790
1791#hc <- hclust(as.dist(Xwer_dist), method = "ward.D")
1792#plot(hc)
1793#
1794# #clust <- cutree(hc, 2)
1795#
1796for(K in 2:30){
1797 #K <- 3
1798 #pamfit <- pam(tdata[-201, ci$selectv], k = K)
1799 pamfit <- pam(as.dist(Xwer_dist), k = K, diss = TRUE)
1800
1801 #table(pamfit$clustering)
1802
1803 SC <- matrix(0, ncol = p, nrow = K)
1804
1805 clustfactor <- pamfit$clustering
1806# for(k in 1:K){
1807# clustk <- which(clustfactor == k)
1808# if(length(clustk) > 0) {
1809# if(length(clustk) > 1) {
1810# SCk <- colSums(synchros09[which(clustfactor == k), ])
1811# } else {
1812# SCk <- synchros09[which(clustfactor == k), ]
1813# }
1814# SC[k, ] <- SC[k, ] + SCk
1815# rm(SCk)
1816# }
1817#}
1818
1819#write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt"))
1820#write.table(clustfactor, file = "~/tmp/clustfactor3.txt")
1821#write.table(clustfactor, file = paste0("~/tmp/clustfactorWER", K, ".txt"))
1822write.table(clustfactor, file = paste0("~/tmp/clustfactor-randomWER", K, ".txt"))
1823}
1824#
1825# # Plots
1826# layout(1)
1827# matplot(t(SC)[48*10 + 1:(48*30), ], type = 'l', ylab = '',col = 1:3, lty = 1)
1828# matplot(t(SC)[48*100 + 1:(48*30), ], type = 'l', ylab = '', col = 1:3, lty = 1)
1829#
1830#
1831#
Electricity Power Curves Clustering