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 BA |
87 | vect <- as.vector(vect) |
88 | matrix(vect[-(1:2)], delta, lscvect) | |
1c6f223e | 89 | } |
1c6f223e | 90 | |
d03c0621 BA |
91 | #fonction smCWT (dans aux.R) |
92 | smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0, nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, | |
93 | lt= 24, dt= 0.5, scalevector ) | |
1c6f223e | 94 | { |
c6556868 BA |
95 | #noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
96 | #scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 | |
d03c0621 BA |
97 | wsp <- Mod(CWT) |
98 | smwsp <- smooth.matrix(wsp, swabs) | |
99 | smsmwsp <- smooth.time(smwsp, tw, dt, scalevector) | |
100 | smsmwsp | |
101 | } | |
1c6f223e | 102 | |
d03c0621 BA |
103 | #dans sowas.R (...donc on ne lisse pas à ce niveau ?) |
104 | smooth.matrix <- function(wt,swabs) | |
105 | { | |
106 | if (swabs != 0) | |
107 | { | |
108 | smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) | |
109 | } else | |
110 | { | |
111 | smwt <- wt | |
112 | } | |
113 | smwt | |
114 | } | |
1c6f223e | 115 | |
d03c0621 BA |
116 | smooth.time <- function(wt,tw,dt,scalevector) |
117 | { | |
118 | smwt <- wt | |
119 | if (tw != 0) | |
120 | { | |
121 | for (i in 1:length(scalevector)) | |
122 | { | |
123 | twi <- as.integer(scalevector[i]*tw/dt) | |
124 | smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) | |
125 | } | |
1c6f223e | 126 | } |
d03c0621 BA |
127 | smwt |
128 | } | |
129 | ||
130 | step2 = function(conso) | |
131 | { | |
132 | #(Benjamin) | |
133 | #à partir de là, "conso" == courbes synchrones | |
134 | n <- nrow(conso) | |
135 | delta <- ncol(conso) | |
136 | ||
137 | #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube] | |
138 | # #NOTE: delta et lscvect pourraient etre gardés à part (communs) | |
139 | ||
140 | #TODO: automatic tune of these parameters ? (for other users) | |
141 | nvoice <- 4 | |
142 | # # noctave4 = 2^13 = 8192 half hours ~ 180 days | |
143 | noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2, tw = 0, noctave = 13) | |
144 | # # 4 here represent 2^5 = 32 half-hours ~ 1 day | |
145 | scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2 | |
146 | lscvect4 <- length(scalevector4) | |
147 | lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect | |
148 | ||
149 | # observations node with CWT | |
150 | Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, scalevector = scalevector4, lt = delta, | |
151 | smooth = FALSE, nvoice = nvoice) | |
3ccd1e39 | 152 | |
d03c0621 BA |
153 | #matrix: |
154 | ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) | |
155 | Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) | |
156 | ||
157 | #NOTE: delta et lscvect pourraient etre gardés à part (communs) | |
158 | for(i in 1:n) | |
159 | Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) | |
160 | #rm(conso, Xcwt4); gc() | |
161 | ||
3ccd1e39 BA |
162 | #Benjamin: FIX is this OK ? |
163 | lscvect = dim(Xcwt4)[2] | |
164 | ||
d03c0621 BA |
165 | ## _.b WER^2 distances ######## |
166 | Xwer_dist <- matrix(0.0, n, n) | |
167 | for(i in 1:(n - 1)) | |
1c6f223e | 168 | { |
3ccd1e39 BA |
169 | #browser() |
170 | ##ERROR là sans FIX lscvect :: delta lscvect --> taille ??! | |
c6556868 BA |
171 | mat1 <- vect2mat(Xcwt2[i,], delta, lscvect) |
172 | ||
173 | for(j in (i + 1):n) | |
d03c0621 | 174 | { |
c6556868 | 175 | mat2 <- vect2mat(Xcwt2[j,], delta, lscvect) |
d03c0621 BA |
176 | num <- Mod(mat1 * Conj(mat2)) |
177 | WX <- Mod(mat1 * Conj(mat1)) | |
178 | WY <- Mod(mat2 * Conj(mat2)) | |
179 | smsmnum <- smCWT(num, scalevector = scalevector4) | |
180 | smsmWX <- smCWT(WX, scalevector = scalevector4) | |
181 | smsmWY <- smCWT(WY, scalevector = scalevector4) | |
182 | wer2 <- sum(colSums(smsmnum)^2) / | |
183 | sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) | |
184 | Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) | |
185 | Xwer_dist[j, i] <- Xwer_dist[i, j] | |
186 | } | |
1c6f223e | 187 | } |
d03c0621 | 188 | diag(Xwer_dist) <- numeric(n) |
c6556868 | 189 | Xwer_dist |
1c6f223e | 190 | } |