temp/WerDist.R

   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
  25 library(Rwave)       # CWT
  26 library(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
 140 cwt.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
 180 wsp <- 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
 280 wco <-
 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
 399 wcs <- 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
 478 rawWSP <- 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
 493 criticalvaluesWSP <- 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
 558 criticalvaluesWCO <- 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
 604 slide <- 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
 663 arealtest <- 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
 690 plotat <- 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
 703 plotcv <- 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
 725 plotcoi <- 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
 736 plotmarks <- 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
 761 plot.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
 767 plotwt <-
 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
 960 createwavesurrogates <- 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
 984 surrwave <- function(x,...)
 985   UseMethod("surrwave")
 986
 987 surrwave.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
1009 surrwave.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
1037 surrwave.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
1071 surrwave.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
1094 invmorlet <- 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
1133 readmatrix <- function(file,M){
1134
1135   A <- matrix(scan(file,comment.char="#"),ncol=M,byrow=TRUE)
1136
1137   A
1138
1139 }
1140
1141
1142 readts <- 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
1164 writematrix <- 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
1175 smooth.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
1186 smooth.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
1204 adjust.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
1218 adjust.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
1229 adjust.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
1243 checkarealsiglevel <- 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
1294 as.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
1338 wtcolors <- 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
1357 createwgn <- function(N,sig,dt){
1358
1359   timeseries<-ts(data=rnorm(N,0,sig),start=0,deltat=dt)
1360
1361   timeseries
1362
1363 }
1364
1365
1366 createar <- 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
1379 rk <- 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
1407 addvalues <- 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
1423 scalematrix <- 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
1437 foldKernel <- 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
1456 kernelBitmap <- 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
1524 kernelRoot <- 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
1549 kernelArea <- 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
1588 tobin <- 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
1599 scaleKernel <- 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
1622 slope <- 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
1706 identifiants <- read.table("identifs.txt")[ ,1]
1707 dates0       <- read.table("datesall.txt")[, 1]
1708 dates        <- as.character(dates0[grep("2009", dates0)])
1709 rm(dates0)
1710
1711 n     <- length(identifiants)
1712 p <- delta <- length(dates)
1713
1714 synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200RC.txt")))
1715 #synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200-random.txt")))
1716
1717 nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing
1718 synchros09[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
1721 imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
1722 synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518])
1723
1724 conso <- synchros09[-201, ];  # series must be on rows
1725 n     <- nrow(conso)
1726 delta <- ncol(conso)
1727
1728 rm(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
1786 load("../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 #
1796 for(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"))
1822 write.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 #