| 1 | #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 |
| 2 | |
| 3 | #(Benjamin) |
| 4 | #à partir de là, "conso" == courbes synchrones |
| 5 | n <- nrow(conso) |
| 6 | delta <- ncol(conso) |
| 7 | |
| 8 | |
| 9 | #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube] |
| 10 | |
| 11 | #TODO: une fonction qui fait lignes 59 à 91 |
| 12 | |
| 13 | #cube: |
| 14 | # Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, |
| 15 | # scalevector = scalevector4, |
| 16 | # lt = delta, smooth = FALSE, |
| 17 | # nvoice = nvoice) # observations node with CWT |
| 18 | # |
| 19 | # #matrix: |
| 20 | # ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) |
| 21 | # #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) |
| 22 | # |
| 23 | # #NOTE: delta et lscvect pourraient etre gardés à part (communs) |
| 24 | # for(i in 1:n) |
| 25 | # Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) |
| 26 | # |
| 27 | # #rm(conso, Xcwt4); gc() |
| 28 | # |
| 29 | # ## _.b WER^2 distances ######## |
| 30 | # Xwer_dist <- matrix(0.0, n, n) |
| 31 | # for(i in 1:(n - 1)){ |
| 32 | # mat1 <- vect2mat(Xcwt2[i,]) |
| 33 | # for(j in (i + 1):n){ |
| 34 | # mat2 <- vect2mat(Xcwt2[j,]) |
| 35 | # num <- Mod(mat1 * Conj(mat2)) |
| 36 | # WX <- Mod(mat1 * Conj(mat1)) |
| 37 | # WY <- Mod(mat2 * Conj(mat2)) |
| 38 | # smsmnum <- smCWT(num, scalevector = scalevector4) |
| 39 | # smsmWX <- smCWT(WX, scalevector = scalevector4) |
| 40 | # smsmWY <- smCWT(WY, scalevector = scalevector4) |
| 41 | # wer2 <- sum(colSums(smsmnum)^2) / |
| 42 | # sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) |
| 43 | # Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) |
| 44 | # Xwer_dist[j, i] <- Xwer_dist[i, j] |
| 45 | # } |
| 46 | # } |
| 47 | # diag(Xwer_dist) <- numeric(n) |
| 48 | # |
| 49 | # save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata") |
| 50 | # save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata") |
| 51 | |
| 52 | |
| 53 | |
| 54 | #lignes 59 à 91 "dépliées" : |
| 55 | Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, |
| 56 | scalevector = scalevector4, |
| 57 | lt = delta, smooth = FALSE, |
| 58 | nvoice = nvoice) # observations node with CWT |
| 59 | |
| 60 | #matrix: |
| 61 | ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) |
| 62 | Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) |
| 63 | |
| 64 | #NOTE: delta et lscvect pourraient etre gardés à part (communs) |
| 65 | for(i in 1:n) |
| 66 | Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) |
| 67 | |
| 68 | #rm(conso, Xcwt4); gc() |
| 69 | |
| 70 | ## _.b WER^2 distances ######## |
| 71 | Xwer_dist <- matrix(0.0, n, n) |
| 72 | for(i in 1:(n - 1)){ |
| 73 | mat1 <- vect2mat(Xcwt2[i,]) |
| 74 | |
| 75 | #NOTE: vect2mat = as.matrix ?! (dans aux.R) |
| 76 | vect2mat <- function(vect){ |
| 77 | vect <- as.vector(vect) |
| 78 | matrix(vect[-(1:2)], delta, lscvect) |
| 79 | } |
| 80 | |
| 81 | for(j in (i + 1):n){ |
| 82 | mat2 <- vect2mat(Xcwt2[j,]) |
| 83 | num <- Mod(mat1 * Conj(mat2)) |
| 84 | WX <- Mod(mat1 * Conj(mat1)) |
| 85 | WY <- Mod(mat2 * Conj(mat2)) |
| 86 | smsmnum <- smCWT(num, scalevector = scalevector4) |
| 87 | smsmWX <- smCWT(WX, scalevector = scalevector4) |
| 88 | smsmWY <- smCWT(WY, scalevector = scalevector4) |
| 89 | wer2 <- sum(colSums(smsmnum)^2) / |
| 90 | sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) |
| 91 | Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) |
| 92 | Xwer_dist[j, i] <- Xwer_dist[i, j] |
| 93 | } |
| 94 | } |
| 95 | diag(Xwer_dist) <- numeric(n) |
| 96 | |
| 97 | #fonction smCWT (dans aux.R) |
| 98 | smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0, |
| 99 | nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, |
| 100 | lt= 24, dt= 0.5, scalevector ) |
| 101 | { |
| 102 | # noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
| 103 | # scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 |
| 104 | wsp <- Mod(CWT) |
| 105 | smwsp <- smooth.matrix(wsp, swabs) |
| 106 | smsmwsp <- smooth.time(smwsp, tw, dt, scalevector) |
| 107 | smsmwsp |
| 108 | } |
| 109 | |
| 110 | #dans sowas.R |
| 111 | smooth.matrix <- function(wt,swabs){ |
| 112 | |
| 113 | if (swabs != 0) |
| 114 | smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) |
| 115 | else |
| 116 | smwt <- wt |
| 117 | |
| 118 | smwt |
| 119 | |
| 120 | } |
| 121 | smooth.time <- function(wt,tw,dt,scalevector){ |
| 122 | |
| 123 | smwt <- wt |
| 124 | |
| 125 | if (tw != 0){ |
| 126 | for (i in 1:length(scalevector)){ |
| 127 | |
| 128 | twi <- as.integer(scalevector[i]*tw/dt) |
| 129 | smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) |
| 130 | |
| 131 | } |
| 132 | } |
| 133 | smwt |
| 134 | } |
| 135 | |
| 136 | #et filter() est dans stats:: |
| 137 | |
| 138 | #cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c |
| 139 | |