| 1 | #include <stdlib.h> |
| 2 | #include <math.h> |
| 3 | #include <stdbool.h> |
| 4 | |
| 5 | #ifndef M_PI |
| 6 | #define M_PI 3.14159265358979323846 |
| 7 | #endif |
| 8 | |
| 9 | // n: number of synchrones, m: length of a synchrone |
| 10 | float computeWerDist(float* s1, float* s2, int n, int m) |
| 11 | { |
| 12 | //TODO: automatic tune of all these parameters ? (for other users) |
| 13 | int nvoice = 4; |
| 14 | //noctave 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones)) |
| 15 | int noctave = 13 |
| 16 | // 4 here represent 2^5 = 32 half-hours ~ 1 day |
| 17 | //NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?) |
| 18 | //R: scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1) |
| 19 | int* scalevector = (int*)malloc( (noctave*nvoice-4 + 1) * sizeof(int)) |
| 20 | for (int i=4; i<=noctave*nvoice; i++) |
| 21 | scalevector[i-4] = pow(2., (float)i/nvoice + 1.); |
| 22 | //condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1 |
| 23 | int s0 = 2; |
| 24 | double w0 = 2*M_PI; |
| 25 | bool scaled = false; |
| 26 | int s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 ) |
| 27 | int totnoct = noctave + as.integer(s0log/nvoice) + 1 |
| 28 | |
| 29 | |
| 30 | |
| 31 | |
| 32 | |
| 33 | ///TODO: continue |
| 34 | |
| 35 | |
| 36 | |
| 37 | computeCWT = function(i) |
| 38 | { |
| 39 | if (verbose) |
| 40 | cat(paste("+++ Compute Rwave::cwt() on serie ",i,"\n", sep="")) |
| 41 | ts <- scale(ts(synchrones[i,]), center=TRUE, scale=scaled) |
| 42 | totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0) |
| 43 | ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)] |
| 44 | #Normalization |
| 45 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) |
| 46 | sqres <- sweep(ts.cwt,2,sqs,'*') |
| 47 | sqres / max(Mod(sqres)) |
| 48 | } |
| 49 | |
| 50 | if (parll) |
| 51 | { |
| 52 | cl = parallel::makeCluster(ncores_clust) |
| 53 | parallel::clusterExport(cl, |
| 54 | varlist=c("synchrones","totnoct","nvoice","w0","s0log","noctave","s0","verbose"), |
| 55 | envir=environment()) |
| 56 | } |
| 57 | |
| 58 | # (normalized) observations node with CWT |
| 59 | Xcwt4 <- |
| 60 | if (parll) |
| 61 | parallel::parLapply(cl, seq_len(n), computeCWT) |
| 62 | else |
| 63 | lapply(seq_len(n), computeCWT) |
| 64 | |
| 65 | if (parll) |
| 66 | parallel::stopCluster(cl) |
| 67 | |
| 68 | Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double") |
| 69 | fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!) |
| 70 | if (verbose) |
| 71 | cat("*** Compute WER distances from CWT\n") |
| 72 | |
| 73 | #TODO: computeDistances(i,j), et répartir les n(n-1)/2 couples d'indices |
| 74 | #là c'est trop déséquilibré |
| 75 | |
| 76 | computeDistancesLineI = function(i) |
| 77 | { |
| 78 | if (verbose) |
| 79 | cat(paste(" Line ",i,"\n", sep="")) |
| 80 | for (j in (i+1):n) |
| 81 | { |
| 82 | #TODO: 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C |
| 83 | num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE) |
| 84 | WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE) |
| 85 | WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE) |
| 86 | wer2 <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) ) |
| 87 | if (parll) |
| 88 | synchronicity::lock(m) |
| 89 | Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2)) |
| 90 | Xwer_dist[j,i] <- Xwer_dist[i,j] |
| 91 | if (parll) |
| 92 | synchronicity::unlock(m) |
| 93 | } |
| 94 | Xwer_dist[i,i] = 0. |
| 95 | } |
| 96 | |
| 97 | parll = (requireNamespace("synchronicity",quietly=TRUE) |
| 98 | && parll && Sys.info()['sysname'] != "Windows") |
| 99 | if (parll) |
| 100 | m <- synchronicity::boost.mutex() |
| 101 | |
| 102 | ignored <- |
| 103 | if (parll) |
| 104 | { |
| 105 | parallel::mclapply(seq_len(n-1), computeDistancesLineI, |
| 106 | mc.cores=ncores_clust, mc.allow.recursive=FALSE) |
| 107 | } |
| 108 | else |
| 109 | lapply(seq_len(n-1), computeDistancesLineI) |
| 110 | Xwer_dist[n,n] = 0. |
| 111 | |
| 112 | mat_dists = matrix(nrow=n, ncol=n) |
| 113 | #TODO: avoid this loop? |
| 114 | for (i in 1:n) |
| 115 | mat_dists[i,] = Xwer_dist[i,] |
| 116 | mat_dists |
| 117 | |