with parallel::export

diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R
index 578b2f3..6090517 100644 (file)
--- a/epclust/R/clustering.R
+++ b/epclust/R/clustering.R
@@ -1,9 +1,9 @@
 # Cluster one full task (nb_curves / ntasks series)
 # Cluster one full task (nb_curves / ntasks series)

-clusteringTask = function(indices_clust)
+clusteringTask = function(indices, ncores)

 {
 {

-       cl_clust = parallel::makeCluster(ncores_clust)
-       parallel::clusterExport(cl_clust,
-               varlist=c("K1","K2","WER"),
+       cl = parallel::makeCluster(ncores)
+       parallel::clusterExport(cl,
+               varlist=c("K1","getCoefs"),

                envir=environment())
        repeat
        {
                envir=environment())
        repeat
        {


@@ -13,44 +13,45 @@ clusteringTask = function(indices_clust)
                                min(nb_series_per_chunk*i,length(indices_clust)), length(indices_clust) )
                        indices_clust[(nb_series_per_chunk*(i-1)+1):upper_bound]
                })
                                min(nb_series_per_chunk*i,length(indices_clust)), length(indices_clust) )
                        indices_clust[(nb_series_per_chunk*(i-1)+1):upper_bound]
                })

-               indices_clust = parallel::parLapply(cl, indices_workers, clusterChunk)
-               # TODO: soft condition between K2 and K1, before applying final WER step
-               if ((WER=="end" && length(indices_clust)==K1) || (WER=="mix" && length(indices_clust)==K2))
+               indices_clust = unlist( parallel::parLapply(cl, indices_workers, function(indices)
+                       computeClusters1(indices, getCoefs, K1)) )
+               if (length(indices_clust) == K1)

                        break
        }
        parallel::stopCluster(cl_clust)
                        break
        }
        parallel::stopCluster(cl_clust)

-       unlist(indices_clust)
+       if (WER == "end")
+               return (indices_clust)
+       #WER=="mix"
+       computeClusters2(indices_clust, K2, getSeries, to_file=TRUE)

 }
 
 }
 

+# Apply the clustering algorithm (PAM) on a coeffs or distances matrix
+computeClusters1 = function(indices, getCoefs, K1)
+       indices[ cluster::pam(getCoefs(indices), K1, diss=FALSE)$id.med ]
+

 # Cluster a chunk of series inside one task (~max nb_series_per_chunk)
 # Cluster a chunk of series inside one task (~max nb_series_per_chunk)

-clusterChunk = function(indices_chunk)
+computeClusters2 = function(indices, K2, getSeries, to_file)

 {
 {

-       coeffs = readCoeffs(indices_chunk)
-       cl = computeClusters(as.matrix(coeffs[,2:ncol(coeffs)]), K1, diss=FALSE)
+       if (is.null(indices))
+       {
+               #get series from file
+       }
+#Puis K-means après WER...

        if (WER=="mix" > 0)
        {
        if (WER=="mix" > 0)
        {

-               curves = computeSynchrones(cl)
+               curves = computeSynchrones(indices)

                dists = computeWerDists(curves)
                dists = computeWerDists(curves)

-               cl = computeClusters(dists, K2, diss=TRUE)
+               indices = computeClusters(dists, K2, diss=TRUE)

        }
        }

-       indices_chunk[cl]
-}
-
-# Apply the clustering algorithm (PAM) on a coeffs or distances matrix
-computeClusters = function(md, K, diss)
-{
-       if (!require(cluster, quietly=TRUE))
-               stop("Unable to load cluster library")
-       cluster::pam(md, K, diss=diss)$id.med
+       if (to_file)
+               #write results to file (JUST series ; no possible ID here)

 }
 
 # Compute the synchrones curves (sum of clusters elements) from a clustering result
 }
 
 # Compute the synchrones curves (sum of clusters elements) from a clustering result

-computeSynchrones = function(indices)
-{
-       colSums( getData(indices) )
-}
+computeSynchrones = function(inds)
+       sapply(seq_along(inds), colMeans(getSeries(inds[[i]]$indices,inds[[i]]$ids)))

 
 

-# Compute the WER distance between the synchrones curves
+# Compute the WER distance between the synchrones curves (in columns)

 computeWerDist = function(curves)
 {
        if (!require("Rwave", quietly=TRUE))
 computeWerDist = function(curves)
 {
        if (!require("Rwave", quietly=TRUE))


@@ -73,7 +74,7 @@ computeWerDist = function(curves)
 
        # (normalized) observations node with CWT
        Xcwt4 <- lapply(seq_len(n), function(i) {
 
        # (normalized) observations node with CWT
        Xcwt4 <- lapply(seq_len(n), function(i) {

-               ts <- scale(ts(curves[i,]), center=TRUE, scale=scaled)
+               ts <- scale(ts(curves[,i]), center=TRUE, scale=scaled)

                totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
                ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
                #Normalization
                totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
                ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
                #Normalization


@@ -88,7 +89,7 @@ computeWerDist = function(curves)
        {
                for (j in (i+1):n)
                {
        {
                for (j in (i+1):n)
                {

-                       #TODO: later, compute CWT here (because not enough storage space for 32M series)
+                       #TODO: later, compute CWT here (because not enough storage space for 200k series)

                        #      'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
                        num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
                        WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
                        #      'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
                        num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
                        WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)

diff --git a/epclust/R/computeCoeffs.R b/epclust/R/computeCoeffs.R
deleted file mode 100644 (file)
index fca3b91..0000000
--- a/epclust/R/computeCoeffs.R
+++ /dev/null
@@ -1,43 +0,0 @@
-computeCoeffs = function(data, index, nb_series_per_chunk, wf)
-{
-       coeffs_chunk = NULL
-       if (is.data.frame(data) && index < nrow(data))
-       {
-               #full data matrix
-               coeffs_chunk = curvesToCoeffs(
-                       data[index:(min(index+nb_series_per_chunk-1,nrow(data))),], wf)
-       }
-       else if (is.function(data))
-       {
-               #custom user function to retrieve next n curves, probably to read from DB
-               coeffs_chunk = curvesToCoeffs( data(rank=(index-1)+seq_len(nb_series_per_chunk)), wf )
-       }
-       else if (exists(data_con))
-       {
-               #incremental connection ; TODO: more efficient way to parse than using a temp file
-               ascii_lines = readLines(data_con, nb_series_per_chunk)
-               if (length(ascii_lines > 0))
-               {
-                       series_chunk_file = ".series_chunk"
-                       writeLines(ascii_lines, series_chunk_file)
-                       coeffs_chunk = curvesToCoeffs( read.csv(series_chunk_file), wf )
-                       unlink(series_chunk_file)
-               }
-       }
-       coeffs_chunk
-}
-
-curvesToCoeffs = function(series, wf)
-{
-       if (!require(wavelets, quietly=TRUE))
-               stop("Couldn't load wavelets library")
-       L = length(series[1,])
-       D = ceiling( log2(L) )
-       nb_sample_points = 2^D
-       #TODO: parallel::parApply() ?!
-       as.data.frame( apply(series, 1, function(x) {
-               interpolated_curve = spline(1:L, x, n=nb_sample_points)$y
-               W = wavelets::dwt(interpolated_curve, filter=wf, D)@W
-               rev( sapply( W, function(v) ( sqrt( sum(v^2) ) ) ) )
-       }) )
-}

diff --git a/epclust/R/main.R b/epclust/R/main.R
index 75041a4..ac4ea8d 100644 (file)
--- a/epclust/R/main.R
+++ b/epclust/R/main.R
@@ -22,6 +22,7 @@
 #' @param ncores_tasks "MPI" number of parallel tasks (1 to disable: sequential tasks)
 #' @param ncores_clust "OpenMP" number of parallel clusterings in one task
 #' @param random Randomize chunks repartition
 #' @param ncores_tasks "MPI" number of parallel tasks (1 to disable: sequential tasks)
 #' @param ncores_clust "OpenMP" number of parallel clusterings in one task
 #' @param random Randomize chunks repartition

+#' @param ... Other arguments to be passed to \code{data} function

 #'
 #' @return A data.frame of the final medoids curves (identifiers + values)
 #'
 #'
 #' @return A data.frame of the final medoids curves (identifiers + values)
 #'


@@ -37,25 +38,40 @@
 #'   + sampleCurves : wavBootstrap de package wmtsa
 #' cl = epclust(getData, K1=200, K2=15, ntasks=1000, nb_series_per_chunk=5000, WER="mix")
 #' @export
 #'   + sampleCurves : wavBootstrap de package wmtsa
 #' cl = epclust(getData, K1=200, K2=15, ntasks=1000, nb_series_per_chunk=5000, WER="mix")
 #' @export

-epclust = function(data, K1, K2, ntasks=1, nb_series_per_chunk=50*K1, min_series_per_chunk=5*K1,
-       wf="haar", WER="end", ncores_tasks=1, ncores_clust=4, random=TRUE)
+epclust = function(series,K1,K2,ntasks=1,nb_series_per_chunk=50*K1,min_series_per_chunk=5*K1,
+       wf="haar",WER="end",ncores_tasks=1,ncores_clust=4,random=TRUE,...)

 {
 {

-       # Check arguments
-       if (!is.data.frame(data) && !is.function(data))
+       # Check/transform arguments
+       bin_dir = "epclust.bin/"
+       dir.create(bin_dir, showWarnings=FALSE, mode="0755")
+       if (!is.function(series))
+       {
+               series_file = paste(bin_dir,"data",sep="")
+               unlink(series_file)
+       }
+       if (is.matrix(series))
+               serialize(series, series_file)
+       else if (!is.function(series))

        {
                tryCatch(
                        {
        {
                tryCatch(
                        {

-                               if (is.character(data))
-                                       data_con = file(data, open="r")
-                               else if (!isOpen(data))
+                               if (is.character(series))
+                                       series_con = file(series, open="r")
+                               else if (!isOpen(series))

                                {
                                {

-                                       open(data)
-                                       data_con = data
+                                       open(series)
+                                       series_con = series

                                }
                                }

+                               serialize(series_con, series_file)
+                               close(series_con)

                        },
                        },

-                       error=function(e) "data should be a data.frame, a function or a valid connection"
+                       error=function(e) "series should be a data.frame, a function or a valid connection"

                )
        }
                )
        }

+       if (!is.function(series))
+               series = function(indices) getDataInFile(indices, series_file)
+       getSeries = series
+

        K1 = toInteger(K1, function(x) x>=2)
        K2 = toInteger(K2, function(x) x>=2)
        ntasks = toInteger(ntasks, function(x) x>=1)
        K1 = toInteger(K1, function(x) x>=2)
        K2 = toInteger(K2, function(x) x>=2)
        ntasks = toInteger(ntasks, function(x) x>=1)


@@ -67,21 +83,22 @@ epclust = function(data, K1, K2, ntasks=1, nb_series_per_chunk=50*K1, min_series
                stop("WER takes values in {'end','mix'}")
 
        # Serialize all wavelets coefficients (+ IDs) onto a file
                stop("WER takes values in {'end','mix'}")
 
        # Serialize all wavelets coefficients (+ IDs) onto a file

-       unlink(".coeffs")
+       coefs_file = paste(bin_dir,"coefs",sep="")
+       unlink(coefs_file)

        index = 1
        nb_curves = 0
        index = 1
        nb_curves = 0

-       nb_coeffs = NA

        repeat
        {
        repeat
        {

-               coeffs_chunk = computeCoeffs(data, index, nb_series_per_chunk, wf)
-               if (is.null(coeffs_chunk))
+               series = getSeries((index-1)+seq_len(nb_series_per_chunk))
+               if (is.null(series))

                        break
                        break

-               writeCoeffs(coeffs_chunk)
+               coeffs_chunk = curvesToCoeffs(series, wf)
+               serialize(coeffs_chunk, coefs_file)

                index = index + nb_series_per_chunk
                nb_curves = nb_curves + nrow(coeffs_chunk)
                index = index + nb_series_per_chunk
                nb_curves = nb_curves + nrow(coeffs_chunk)

-               if (is.na(nb_coeffs))
-                       nb_coeffs = ncol(coeffs_chunk)-1

        }
        }

+       getCoefs = function(indices) getDataInFile(indices, coefs_file)
+######TODO: if DB, array rank --> ID at first retrieval, when computing coeffs; so:: NO use of IDs !

 
        if (nb_curves < min_series_per_chunk)
                stop("Not enough data: less rows than min_series_per_chunk!")
 
        if (nb_curves < min_series_per_chunk)
                stop("Not enough data: less rows than min_series_per_chunk!")


@@ -95,16 +112,17 @@ epclust = function(data, K1, K2, ntasks=1, nb_series_per_chunk=50*K1, min_series
                upper_bound = ifelse( i<ntasks, min(nb_series_per_task*i,nb_curves), nb_curves )
                indices[((i-1)*nb_series_per_task+1):upper_bound]
        })
                upper_bound = ifelse( i<ntasks, min(nb_series_per_task*i,nb_curves), nb_curves )
                indices[((i-1)*nb_series_per_task+1):upper_bound]
        })

-       library(parallel, quietly=TRUE)

        cl_tasks = parallel::makeCluster(ncores_tasks)
        parallel::clusterExport(cl_tasks,
        cl_tasks = parallel::makeCluster(ncores_tasks)
        parallel::clusterExport(cl_tasks,

-               varlist=c("K1","K2","WER","nb_series_per_chunk","ncores_clust"),#TODO: pass also
-                                               #nb_coeffs...and filename (in a list... ?)
+               varlist=c("getSeries","getCoefs","K1","K2","WER","nb_series_per_chunk","ncores_clust"),

                envir=environment())
                envir=environment())

+       #1000*K1 (or K2) indices (or NOTHING--> series on file)

        indices = parallel::parLapply(cl_tasks, indices_tasks, clusteringTask)
        parallel::stopCluster(cl_tasks)
 
        indices = parallel::parLapply(cl_tasks, indices_tasks, clusteringTask)
        parallel::stopCluster(cl_tasks)
 

-       # Run step1+2 step on resulting ranks
-       indices = clusterChunk(indices, K1, K2)
-       return (list("indices"=indices, "medoids"=getSeries(data, indices)))
+       #Now series must be retrieved from synchrones_file, and have no ID
+       getSeries = function(indices, ids) getDataInFile(indices, synchrones_file)
+
+       # Run step2 on resulting indices or series (from file)
+       computeClusters2(indices=if (WER=="end") indices else NULL, K2, to_file=FALSE)

 }
 }

diff --git a/epclust/R/utils.R b/epclust/R/utils.R
index e0f25ec..7083674 100644 (file)
--- a/epclust/R/utils.R
+++ b/epclust/R/utils.R
@@ -31,3 +31,15 @@ getSeries(data, rank=NULL, id=NULL)
 {
        #TODO:
 }
 {
        #TODO:
 }

+
+curvesToCoeffs = function(series, wf)
+{
+       L = length(series[1,])
+       D = ceiling( log2(L) )
+       nb_sample_points = 2^D
+       apply(series, 1, function(x) {
+               interpolated_curve = spline(1:L, x, n=nb_sample_points)$y
+               W = wavelets::dwt(interpolated_curve, filter=wf, D)@W
+               rev( sapply( W, function(v) ( sqrt( sum(v^2) ) ) ) )
+       })
+}