Commit | Line | Data |
---|---|---|
d03c0621 BA |
1 | #point avec Jairo: |
2 | #rentrer dans code C cwt continue Rwave | |
3 | #passer partie sowas à C | |
4 | #fct qui pour deux series (ID, medoides) renvoie distance WER (Rwave ou à moi) | |
5 | #transformee croisee , smoothing lissage 3 composantes , + calcul pour WER | |
6 | #attention : code fait pour des series temps desynchronisees ! (deltat, dt == 1,2 ...) | |
7 | #determiner nvoice noctave (entre octave + petit et + grand) | |
8 | ||
1c6f223e BA |
9 | library("Rwave") |
10 | ||
dc1aa85a | 11 | #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 |
1c6f223e | 12 | #TODO: bout de code qui calcule les courbes synchrones après étapes 1+2 à partir des ID médoïdes |
dc1aa85a | 13 | |
1c6f223e BA |
14 | #toCWT: (aux) |
15 | ##NOTE: renvoie une matrice 3D | |
c6556868 BA |
16 | toCWT <- function(X, sw=0, tw=0, swabs=0, nvoice=12, noctave=5, s0=2, w0=2*pi, |
17 | lt=24, dt=0.5, spectra=FALSE, smooth=TRUE, scaled=FALSE, scalevector) | |
d03c0621 BA |
18 | { |
19 | noctave <- adjust.noctave(lt, dt, s0, tw, noctave) | |
20 | if(missing(scalevector)) | |
21 | scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 | |
22 | res <- lapply(1:nrow(X), function(n) { | |
23 | tsX <- ts( X[n,] ) | |
24 | tsCent <- tsX - mean(tsX) | |
25 | if(scaled) | |
26 | tsCent <- ts(scale(tsCent)) | |
27 | tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0) | |
28 | tsCent.cwt | |
29 | }) | |
30 | if( spectra ) | |
31 | res <- lapply(res, function(l) Mod(l)^2 ) | |
32 | if( smooth ) | |
33 | res <- lapply(res, smCWT, swabs = swabs, tw = tw, dt = dt, scalevector = scalevector) | |
34 | resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]), length(res))) | |
35 | for( l in 1:length(res) ) | |
36 | resArray[ , , l] <- res[[l]] | |
37 | resArray | |
1c6f223e BA |
38 | } |
39 | ||
c6556868 BA |
40 | #from sowas |
41 | adjust.noctave <- function(N,dt,s0,tw,noctave) | |
42 | { | |
43 | if (tw>0) | |
44 | { | |
45 | dumno <- as.integer((log(N*dt)-log(2*tw*s0))/log(2)) | |
46 | if (dumno<noctave) | |
47 | { | |
48 | cat("# noctave adjusted to time smoothing window \n") | |
49 | noctave <- dumno | |
50 | } | |
51 | } | |
52 | noctave | |
53 | } | |
54 | ||
d03c0621 BA |
55 | #from sowas |
56 | cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi) | |
57 | { | |
58 | if (class(ts)!="ts") | |
59 | stop("# 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") | |
60 | ||
61 | t=time(ts) | |
62 | dt=t[2]-t[1] | |
63 | s0unit=s0/dt*w0/(2*pi) | |
64 | s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) | |
65 | if (s0log<1) | |
66 | { | |
67 | cat(paste("# s0unit = ",s0unit,"\n",sep="")) | |
68 | cat(paste("# s0log = ",s0log,"\n",sep="")) | |
69 | cat("# s0 too small for w0! \n") | |
70 | } | |
71 | totnoct=noctave+as.integer(s0log/nvoice)+1 | |
dc1aa85a | 72 | |
d03c0621 BA |
73 | #cwt from package Rwave |
74 | totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0) | |
75 | ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)] | |
dc1aa85a | 76 | |
d03c0621 BA |
77 | #Normalization |
78 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) | |
79 | smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) | |
dc1aa85a | 80 | |
d03c0621 | 81 | ts.cwt*smat |
dc1aa85a BA |
82 | } |
83 | ||
d03c0621 | 84 | #NOTE: vect2mat = as.matrix ?! (dans aux.R) |
c6556868 | 85 | vect2mat <- function(vect, delta, lscvect) |
1c6f223e | 86 | { |
d03c0621 | 87 | vect <- as.vector(vect) |
c6556868 BA |
88 | |
89 | print(delta) | |
90 | print(lscvect) | |
91 | print(delta * lscvect) | |
92 | browser() | |
93 | ||
94 | ||
d03c0621 | 95 | matrix(vect[-(1:2)], delta, lscvect) |
1c6f223e | 96 | } |
1c6f223e | 97 | |
d03c0621 BA |
98 | #fonction smCWT (dans aux.R) |
99 | smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0, nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, | |
100 | lt= 24, dt= 0.5, scalevector ) | |
1c6f223e | 101 | { |
c6556868 BA |
102 | #noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
103 | #scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 | |
d03c0621 BA |
104 | wsp <- Mod(CWT) |
105 | smwsp <- smooth.matrix(wsp, swabs) | |
106 | smsmwsp <- smooth.time(smwsp, tw, dt, scalevector) | |
107 | smsmwsp | |
108 | } | |
1c6f223e | 109 | |
d03c0621 BA |
110 | #dans sowas.R (...donc on ne lisse pas à ce niveau ?) |
111 | smooth.matrix <- function(wt,swabs) | |
112 | { | |
113 | if (swabs != 0) | |
114 | { | |
115 | smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) | |
116 | } else | |
117 | { | |
118 | smwt <- wt | |
119 | } | |
120 | smwt | |
121 | } | |
1c6f223e | 122 | |
d03c0621 BA |
123 | smooth.time <- function(wt,tw,dt,scalevector) |
124 | { | |
125 | smwt <- wt | |
126 | if (tw != 0) | |
127 | { | |
128 | for (i in 1:length(scalevector)) | |
129 | { | |
130 | twi <- as.integer(scalevector[i]*tw/dt) | |
131 | smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) | |
132 | } | |
1c6f223e | 133 | } |
d03c0621 BA |
134 | smwt |
135 | } | |
136 | ||
137 | step2 = function(conso) | |
138 | { | |
139 | #(Benjamin) | |
140 | #à partir de là, "conso" == courbes synchrones | |
141 | n <- nrow(conso) | |
142 | delta <- ncol(conso) | |
143 | ||
144 | #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube] | |
145 | # #NOTE: delta et lscvect pourraient etre gardés à part (communs) | |
146 | ||
147 | #TODO: automatic tune of these parameters ? (for other users) | |
148 | nvoice <- 4 | |
149 | # # noctave4 = 2^13 = 8192 half hours ~ 180 days | |
150 | noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2, tw = 0, noctave = 13) | |
151 | # # 4 here represent 2^5 = 32 half-hours ~ 1 day | |
152 | scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2 | |
153 | lscvect4 <- length(scalevector4) | |
154 | lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect | |
155 | ||
156 | # observations node with CWT | |
157 | Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, scalevector = scalevector4, lt = delta, | |
158 | smooth = FALSE, nvoice = nvoice) | |
c6556868 | 159 | browser() |
d03c0621 BA |
160 | #matrix: |
161 | ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) | |
162 | Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) | |
163 | ||
164 | #NOTE: delta et lscvect pourraient etre gardés à part (communs) | |
165 | for(i in 1:n) | |
166 | Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) | |
167 | #rm(conso, Xcwt4); gc() | |
168 | ||
169 | ## _.b WER^2 distances ######## | |
170 | Xwer_dist <- matrix(0.0, n, n) | |
171 | for(i in 1:(n - 1)) | |
1c6f223e | 172 | { |
d03c0621 | 173 | |
c6556868 BA |
174 | |
175 | ||
176 | ||
177 | ||
178 | ||
179 | ##ERROR là :: delta lscvect --> taille ??! | |
180 | mat1 <- vect2mat(Xcwt2[i,], delta, lscvect) | |
181 | ||
182 | for(j in (i + 1):n) | |
d03c0621 | 183 | { |
c6556868 | 184 | mat2 <- vect2mat(Xcwt2[j,], delta, lscvect) |
d03c0621 BA |
185 | num <- Mod(mat1 * Conj(mat2)) |
186 | WX <- Mod(mat1 * Conj(mat1)) | |
187 | WY <- Mod(mat2 * Conj(mat2)) | |
188 | smsmnum <- smCWT(num, scalevector = scalevector4) | |
189 | smsmWX <- smCWT(WX, scalevector = scalevector4) | |
190 | smsmWY <- smCWT(WY, scalevector = scalevector4) | |
191 | wer2 <- sum(colSums(smsmnum)^2) / | |
192 | sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) | |
193 | Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) | |
194 | Xwer_dist[j, i] <- Xwer_dist[i, j] | |
195 | } | |
1c6f223e | 196 | } |
d03c0621 | 197 | diag(Xwer_dist) <- numeric(n) |
c6556868 | 198 | Xwer_dist |
1c6f223e | 199 | } |