first commit

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
index 0000000..37650a7
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,8 @@
+.Rhistory
+.RData
+*.o
+*.so
+*.dll
+
+/*.zip
+/SYNCLUST_init/*

diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100755 (executable)
index 0000000..b9abed1
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,15 @@
+Package: synclust
+Type: Package
+Version: 0.1.1
+Date: 2013-01-31
+Title: Delimiting synchronous population areas
+Author: Benjamin Auder, Christophe Giraud
+Maintainer: Benjamin Auder <Benjamin.Auder@gmail.com>
+Depends: R (>= 2.14.1), mvtnorm
+Suggests: kernlab
+Description: Provide two methods to cluster species by regions,
+        using temporal variations and/or geographic coordinates. 
+        The resulting areas (should) have synchronous variations.
+License: GPL (>= 3)
+LazyData: yes
+LazyLoad: yes

diff --git a/NAMESPACE b/NAMESPACE
new file mode 100755 (executable)
index 0000000..729d29a
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,8 @@
+# Export all user-level R functions
+export (findSyncVarRegions, drawMapWithSites, 
+               drawNeighborhoodGraph, plotCurves, .Last.lib)
+
+# Import all packages listed as Imports or Depends
+#import (methods)
+
+useDynLib(synclust)

diff --git a/R/clustering.R b/R/clustering.R
new file mode 100644 (file)
index 0000000..e5bdbf1
--- /dev/null
+++ b/R/clustering.R
@@ -0,0 +1,35 @@
+#main function (choice between kmeans and hierarchical clustering)
+getClusters = function(distances, method, K)
+{
+       clusts = c()
+       if (method=="KM")
+       {
+               nstart = 10 #number of kmeans random restarts
+               maxiter = 100 #maximum iterations count in each km run
+               clusts = .Call("kmeansWithDistances", distances, K, nstart, maxiter)
+       }
+       else if (method=="HC")
+       {
+               #simple hierarchical clustering using ECT distances
+               hct = hclust(as.dist(distances),method="ward.D")
+               clusts = cutree(hct, K)
+       }
+       return (clusts)
+}
+
+# renumbering step (post-processing after clustering)
+reordering = function(clusts)
+{
+       newCl = clusts
+       maxInd = max(clusts)
+       counter = 1
+       for (i in 1:maxInd)
+       {
+               if (sum(clusts == i) > 0)
+               {
+                       newCl[clusts == i] = counter
+                       counter = counter + 1
+               }
+       }
+       return (newCl)
+}

diff --git a/R/distances.R b/R/distances.R
new file mode 100644 (file)
index 0000000..8063cec
--- /dev/null
+++ b/R/distances.R
@@ -0,0 +1,80 @@
+#build similarity matrix W (NOTE : sparse matrix ==> optimizations later)
+getSimilarityMatrix = function(NI)
+{
+       # using a local sigma would be nice, but break W symmetry,
+       # which cannot easily be repaired then (??!)
+       # ==> we use a global sigma, with a very simple heuristic
+       
+       n = length(NI$ix)
+       distances = c()
+       for (i in 1:n) distances = c(distances,NI$ds[[i]])
+       distances = unique(distances)
+       sigma2 = median(distances)^2 #for example...
+
+       W = matrix(0.0,nrow=n,ncol=n)
+       for (i in 1:n)
+               W[ i, NI$ix[[i]] ] = exp( - NI$ds[[i]]^2 / sigma2 )
+
+       return (W)
+}
+
+#epsilon constant, used as a zero threshold
+EPS = 100 * .Machine$double.eps
+
+#Moore-Penrose pseudo inverse
+mppsinv = function(M)
+{
+       s = svd(M)
+       sdiag = s$d ; sdiag[sdiag < EPS] = Inf
+       p = min(nrow(M),ncol(M))
+       sdiag = diag(1.0 / sdiag, p)
+       return ((s$v) %*% sdiag %*% t(s$u))
+}
+
+#get distance matrix from data and similarity : Commute Time
+getECTDistances = function(NI)
+{
+       n = length(NI$ix) ; seqVect = 1:n
+       if (n <= 1) return (0.0) #nothing to do...
+       
+       #get laplacian (...inverse) :
+       W = getSimilarityMatrix(NI)
+       invLap = mppsinv(diag(rowSums(W)) - W)
+
+       #...and distances
+       ectd = matrix(0.0, nrow=n, ncol=n)
+       for (ij in 1:n)
+       {
+               ectd[ij,] = ectd[ij,] + invLap[ij,ij]
+               ectd[,ij] = ectd[,ij] + invLap[ij,ij]
+       }
+       ectd = ectd - 2*invLap
+       return (ectd)
+}
+
+# Call Dijsktra algorithm on every vertex to build distances matrix
+getShortestPathDistances = function(NI)
+{
+       n = length(NI$ix)
+       distancesIn = matrix(NA,nrow=n,ncol=n)
+       for (i in 1:n)
+               distancesIn[i,NI$ix[[i]]] = NI$ds[[i]]
+
+       distancesOut = matrix(nrow=n, ncol=n)
+       for (i in 1:n)
+               distancesOut[i,] = .Call("dijkstra", distancesIn, i)
+       return (distancesOut)
+}
+
+## MAIN CALL to get distances matrix
+getDistances = function(dtype, NI)
+{
+       distances = matrix()
+       if (dtype=="spath")
+               distances = getShortestPathDistances(NI)
+       else if (dtype=="ectd")
+               distances = getECTDistances(NI)
+       
+       diag(distances) = 0.0 #distances to self are zero
+       return (distances)
+}

diff --git a/R/graphics.R b/R/graphics.R
new file mode 100644 (file)
index 0000000..9ce8a3a
--- /dev/null
+++ b/R/graphics.R
@@ -0,0 +1,40 @@
+#draw (France or...) map with all sites of colors 'cols'
+drawMapWithSites = function(M, cols=rep(1,nrow(M)))
+{
+       xMin = range(M[,1])[1]
+       xMax = range(M[,1])[2]
+       yMin = range(M[,2])[1]
+       yMax = range(M[,2])[2]
+       par(mar=c(2,2,2,2))
+       plot(0,0,xlim=c(xMin,xMax),ylim=c(yMin,yMax),col="white")
+       #plot by color groups (limited to integers)
+       maxColor = max(cols)
+       for (i in 1:maxColor)
+       {
+               indices = (1:nrow(M))[cols==i]
+               if (length(indices) > 0)
+                       points(M[indices,1],M[indices,2],col=i,xtitle=NULL)
+       }
+}
+
+#draw neighborhoods graph on top of a country map (or any other map)
+drawNeighborhoodGraph = function(M, NI)
+{
+       for (i in 1:length(NI))
+       {
+               for (j in NI[[i]])
+                       lines(c(M[i,1],M[j,1]),c(M[i,2],M[j,2]))
+       }
+}
+
+#plot a matrix of curves (in rows)
+plotCurves = function(M, cols=rep(1,nrow(M)))
+{
+       n = nrow(M)
+       rg = c(min(M),max(M)) #range for plotting
+       for (i in 1:n)
+       {
+               plot(M[i,],col=cols[i],ylim=rg,type="l")
+               if (i < n) par(new=TRUE)
+       }
+}

diff --git a/R/main.R b/R/main.R
new file mode 100644 (file)
index 0000000..0160b58
--- /dev/null
+++ b/R/main.R
@@ -0,0 +1,207 @@
+#example of "not too bad" parameters
+#~ k=10
+#~ alpha=0.1 
+#~ gmode=1 
+#~ K = 5 
+#~ dtype = "spath"
+#~ cmeth = "HC"
+#~ pcoef=??
+#~ h=??
+#~ eps=??
+#~ maxit=??
+
+#MAIN FUNCTION : direct clustering from a neighborhoods graph, or get regions
+#from (Poisson) distribution parameters optimization, using convex relaxation.
+findSyncVarRegions = function(
+       method, #global method: "direct" or "convex"
+       M, #matrix of observations in rows, the two last columns 
+          #corresponding to geographic coordinates; 
+          #set to NULL to use our initial dataset (625 rows / 9 years)
+       k, #number of neighbors
+       alpha, #weight parameter for intra-neighborhoods distance computations
+              #0 = take only geographic coordinates into account
+              #1 = take only observations over the years into account
+              #in-between : several levels of compromise
+              #-1 or any negative value : use a heuristic to choose alpha
+       gmode, #0 = reduced [mutual] kNN; 1 = augmented kNN; (symmetric)
+              #2 = normal kNN; 3 = one NN in each quadrant; (NON-symmetric)
+                  #NOTE: gmode==3 automatically sets k==4 (at most!)
+       K, #number of clusters
+       dtype, #distance type, in {"simple","spath","ectd"}.
+              #NOTE: better avoid "simple" if gmode>=2
+       cmeth, #clustering method, in {"KM","HC","spec"} for k-means (distances based) 
+              #or hierarchical clustering, or spectral clustering (only if gmode>=2)
+       pcoef=1.0, #penalty value for convex optimization
+       h=1e-3, #step in the min LL algorithm
+       eps=1e-3, #threshold to stop min.LL iterations
+       maxit=1e3, #maximum number of iterations in the min LL algo
+       showLL=TRUE, #print trace of log-likelihood evolution
+       disp=TRUE #true for interactive display (otherwise nothing gets plotted)
+) {
+       #get matrix M if not directly provided
+       if (is.null(M))
+       {
+               data("example", package="synclust")
+               M = synclust_sample
+       }
+       if (is.character(M))
+               M = as.matrix(read.table(M))
+
+       n = nrow(M)
+       m = ncol(M)
+
+       #pretreatment for neighborhoods search: standardize M columns
+       #TODO: maybe apply only on coordinates columns ?
+       std = standardize(M)
+
+       #get neighborhoods [FALSE because NOT simpleDists; see C code]
+       NI = .Call("getNeighbors", std$M, k, alpha, gmode, FALSE)
+
+       #optional intermediate display : map + graph (monocolor)
+       if (disp)
+               promptForMapDisplay("interm", M[,(m-1):m], NIix=NI$ix)
+
+       clusters = rep(1,n)
+       distances = matrix(NA,nrow=n,ncol=n)
+       cxpar = list()
+
+       ## DIRECT CLUSTERING ##
+       if (method=="direct")
+       {
+               if (gmode >= 2)
+                       stop("'gmode' must be 0 or 1 for direct clustering")
+               if (dtype=="simple")
+                       stop("'dtype' cannot be set to \"simple\" for direct (graph!) clustering")
+
+               #find connected components in the graph defined by NI
+               cc = reordering(.Call("getConnectedComponents", NI$ix))
+               nbC = max(cc)
+               if (nbC > 10)
+                       stop(paste("ABORT: too many connex components (found ",nbC,")",sep=''))
+               if (nbC > 1)
+                       print(paste("*** WARNING:",nbC,"connex components ***"))
+               clusters = cc
+
+               #for each connected component...
+               for (i in 1:nbC)
+               {
+                       indices = (1:n)[cc == i]
+                       nc = length(indices)
+                       if (nc <= 1)
+                               next #nothing to do with such a small component
+                       
+                       if (nbC > 1)
+                       {
+                               doClust = readline(paste(">>> cluster current component of cardinal",nc,"/",n,"? (y/n)\n"))
+                               if (doClust == "y")
+                                       K = readline(">>> into how many groups ? (int >= 2)\n")
+                               else
+                                       next
+                       }
+                       
+                       #0] remap NI in current connex component
+                       locNI = remapNeighbors(NI, indices)
+                       
+                       #1] determine similarity and distance matrices (e.g. using a random walk)
+                       locDists = getDistances(dtype, locNI)
+                       distances[indices,indices] = locDists
+                       
+                       #2] cluster data inside connex component according to distance matrix
+                       locClusters = getClusters(locDists, cmeth, K)
+                       maxInd = max(clusters)
+                       clusters[indices] = locClusters + maxInd #avoid indices overlaps
+               }
+       }
+
+       ## CONVEX RELAXATION ##
+       else if (method=="convex")
+       {               
+               #preliminary: remove NA's by averaging over each serie's values
+               M = replaceNAs(M)
+               
+               #use NI$ix and matrix M to estimate initial parameters,
+               #and then iterate until convergence to get f + theta
+               #theta == mean observations count at each site s
+               #f == estimated variations at each site ("time-series" of T points)
+               cxpar = .Call("getVarsWithConvexOptim", 
+                       M[,1:(m-2)], NI$ix, pcoef, h, eps, maxit, (gmode <= 1), showLL)
+               f = cxpar$f #the only one we use (others can be checked by user later)
+               
+               #cluster "time-series" f, using simple kmeans/HC, spect.clust, 
+               #or [in a graph] KM or HC, after redefining a NI (using f only)
+               
+               if (dtype=="simple")
+               {
+                       #use R core functions
+                       if (cmeth=="KM")
+                               clusters = kmeans(f, K, iter.max=100, nstart=10)$cluster
+                       else if (cmeth=="HC")
+                       {
+                               hct = hclust(dist(f), method="ward")
+                               clusters = cutree(hct, K)
+                       }
+                       else if (cmeth=="spec")
+                       {
+                               require(kernlab)
+                               clusters = as.integer(specc(f, K, kpar="automatic"))
+                       }
+               }
+               
+               else
+               {
+                       # recompute NI from repaired/smoothed data [simpleDists=TRUE, no graph dists]
+                       #NOTE: gmode==1, augmented kNN (arbitrary, but should be symmetric)
+                       NI = .Call("getNeighbors", f, k, alpha, 1, TRUE)
+                       
+                       #find connected components in the graph defined by NI
+                       cc = reordering(.Call("getConnectedComponents", NI$ix))
+                       
+                       nbC = max(cc)
+                       if (nbC > 10) 
+                               stop(paste("ABORT: too many connex components (found ",nbC,")",sep=''))
+                       if (nbC > 1) 
+                               print(paste("*** WARNING:",nbC,"connex components ***"))
+                       clusters = cc
+                       
+                       #for each connected component...
+                       for (i in 1:nbC)
+                       {
+                               indices = (1:n)[cc == i]
+                               nc = length(indices)
+                               if (nc <= 1)
+                                       next #nothing to do with such a small component
+                               
+                               if (nbC > 1)
+                               {
+                                       doClust = readline(paste(">>> cluster current component of cardinal",nc,"/",n,"? (y/n)\n"))
+                                       if (doClust == "y")
+                                               K = readline(">>> into how many groups ? (int >= 2)\n")
+                                       else
+                                               next
+                               }
+                               
+                               #0] remap NI in current connex component
+                               locNI = remapNeighbors(NI, indices)
+                               
+                               #1] determine similarity and distance matrices (e.g. using a random walk)
+                               locDists = getDistances(dtype, locNI)
+                               distances[indices,indices] = locDists
+                               
+                               #2] cluster data inside connex component according to distance matrix
+                               locClusters = getClusters(locDists, cmeth, K)
+                               maxInd = max(clusters)
+                               clusters[indices] = locClusters + maxInd #avoid indices overlaps
+                       }
+               }
+       }
+
+       clusters = reordering(clusters)
+       #optional final display : map with clusters colors
+       if (disp)
+               promptForMapDisplay("final", M[,(m-1):m], clusters=clusters)
+
+       #give back matrix M as given to the function
+       M = destandardize(std)
+
+       return (list("M"=M, "NI"=NI, "dists"=distances, "clusts"=clusters, "cxpar"=cxpar))
+}

diff --git a/R/main.utils.R b/R/main.utils.R
new file mode 100644 (file)
index 0000000..6d25c08
--- /dev/null
+++ b/R/main.utils.R
@@ -0,0 +1,85 @@
+#various util functions for the main program
+
+#preliminary: replace NA's by averaging over each serie's values
+#TODO: find a better way to handle missing values
+replaceNAs = function(M)
+{
+       n = nrow(M)
+       m = ncol(M)
+       res = M
+       for (i in 1:n)
+       {
+               avg = mean(M[i,1:(m-2)] [!is.na(M[i,1:(m-2)])])
+               res[i,1:(m-2)] [is.na(M[i,1:(m-2)])] = avg
+       }
+       return (res)
+}
+
+#standardize matrix M (remove mean, divide by standard deviation)
+standardize = function(M)
+{
+       avgM = colMeans(M, na.rm = TRUE)
+       stdevs = sqrt( unlist( apply(M, 2, var, na.rm=TRUE) ) )
+       res = t(M) - avgM
+       res = t(res / stdevs)
+       return (list("M"=res,"avg"=avgM,"stv"=stdevs))
+}
+
+#opposite of the previous function: get back M from standardized form
+destandardize = function(std)
+{
+       M = std$M
+       M = t(M) * std$stv
+       M = t(M + std$avg)
+       return (M)
+}
+
+#remap neighbors into some connex component
+remapNeighbors = function(NI, indices)
+{
+       revIndices = rep(NA, length(NI))
+       nc = length(indices)
+       for (ii in 1:nc)
+               revIndices[ indices[ii] ] = ii
+       locNI = list("ix"=as.list(rep(NA,nc)), "ds"=as.list(rep(NA,nc)))
+       for (ii in 1:nc)
+       {
+               locNI$ix[[ii]] = revIndices[ NI$ix[[ indices[ii] ]] ]
+               locNI$ds[[ii]] = NI$ds[[ indices[ii] ]]
+       }
+       return (locNI)
+}
+
+#check graph connexity
+getConnectedComponents = function(NIix)
+{
+       return (.Call("getConnectedComponents", NIix));
+}
+
+#auxiliary function to display clustering information
+promptForMapDisplay = function(stage, coordsM, NIix=NULL, clusters=NULL)
+{
+       if (is.null(clusters))
+               clusters = rep(1, nrow(coordsM))
+
+       shouldDisplay = ""
+       if (stage == "interm")
+               shouldDisplay = readline(">>> show intermediate map of neighborhoods ? (y/n)\n")
+       else if (stage == "final")
+       {
+               shouldDisplay = readline(
+">>> show final map of clusters ? (y/n) \
+NOTE: can be plotted later, see '? drawMapWithSites'\n")
+       }
+
+       if (shouldDisplay == "y")
+       {
+               drawMapWithSites(coordsM, clusters)
+               if (!is.null(NIix))
+                       drawNeighborhoodGraph(coordsM,NIix)
+               print("Please press 'enter' to continue")
+               readline()
+               if (!is.null(dev.list()))
+                       dev.off()
+       }
+}

diff --git a/R/tests/helpers.R b/R/tests/helpers.R
new file mode 100644 (file)
index 0000000..0650394
--- /dev/null
+++ b/R/tests/helpers.R
@@ -0,0 +1,14 @@
+checkEquals = function(target, current, 
+       tolerance=.Machine$double.eps^0.5, msg="", ...)
+{
+       #all.equal.numeric ?
+       result = all.equal(target, current, tolerance=tolerance, ...)
+       if (result != TRUE)
+               stop(msg)
+}
+
+checkTrue = function(b, msg="")
+{
+       if (!b)
+               stop(msg)
+}

diff --git a/R/tests/runAll.R b/R/tests/runAll.R
new file mode 100644 (file)
index 0000000..610168d
--- /dev/null
+++ b/R/tests/runAll.R
@@ -0,0 +1,20 @@
+source("helpers.R")
+source("t.clustering.R")
+source("t.connexity.R")
+source("t.utils.R")
+
+dyn.load("../../src/synclust.so")
+
+functions = c(lsf.str())
+for (func in functions)
+{
+       #ou test length(grep("test.", function)) > 0
+       if (nchar(func) > 5 && substr(func, 1, 5) == "test.")
+       {
+               print(paste("run",func))
+               eval(call(func))
+       }
+}
+
+#sample call for full package :
+#t = findSyncVarRegions(method="convex",M=NULL,k=10,alpha=0.0,gmode=1,K=5,dtype="spath",cmeth="HC",pcoef=2.2,h=5e-4,eps=1e-5,maxit=3e3,showLL=TRUE,disp=TRUE)

diff --git a/R/tests/t.clustering.R b/R/tests/t.clustering.R
new file mode 100644 (file)
index 0000000..013ddf2
--- /dev/null
+++ b/R/tests/t.clustering.R
@@ -0,0 +1,105 @@
+#test several clustering methods on iris dataset (setosa should be found)
+test.clustering1 = function()
+{
+       data(iris)
+
+       #get neighborhoods from data [25 is high, but shouldn't be lower to have 1 connex comp.]
+       NI = .Call("getNeighbors", as.matrix(iris[,1:4]), 25, 0.0, 1, TRUE)
+
+       for (dtype in c("spath"))#,"ectd")) #bug: TODO
+       {
+               #get distances from neighborhoods; should be OK for all distances 
+               #except "simple" (which is treated as a special case with built-in R funcs)
+               distances = synclust:::getDistances(dtype, NI)
+
+               for (cmeth in c("KM","HC"))
+               {
+                       #finally, get clusters
+                       clusters = synclust:::getClusters(distances, cmeth, K=3)
+                       #check that cluster 'setosa' is pure and separated
+                       uqclust = unique(clusters[1:50])
+                       checkTrue(length(uqclust) == 1)
+                       checkTrue(! uqclust[1] %in% clusters[51:150])
+               }
+       }
+}
+
+#test several parameters agencements on custom non-isotropic gaussian dataset (2D)
+test.clustering2 = function()
+{
+       clustSize = 33
+       
+       require(mvtnorm)
+       set.seed(32)
+       gaussian1 = rmvnorm(clustSize, mean=c(-4.0,-6.0), sigma=matrix(c(1.0,0.7,0.7,1.0),nrow=2))
+       gaussian2 = rmvnorm(clustSize, mean=c(0.0,0.0), sigma=matrix(c(1.0,0.0,0.0,1.0),nrow=2))
+       gaussian3 = rmvnorm(clustSize, mean=c(4.0,-6.0), sigma=matrix(c(1.0,-0.7,-0.7,1.0),nrow=2))
+       M = rbind(gaussian1, rbind(gaussian2, gaussian3))
+       
+       #get neighborhoods from data [25 is high, but shouldn't be much lower to have 1 connex comp.]
+       NI = .Call("getNeighbors", M, 25, 0.0, 1, TRUE)
+       
+       for (dtype in c("spath"))#,"ectd")) #TODO
+       {
+               #get distances from neighborhoods; should be OK for all distances 
+               #except "simple" (which is treated as a special case with built-in R funcs)
+               distances = synclust:::getDistances(dtype, NI)
+               
+               for (cmeth in c("KM","HC"))
+               {
+                       #finally, get clusters
+                       clusters = synclust:::getClusters(distances, cmeth, K=3)
+                       
+                       #soft check, because have to succeed at each run
+                       srt = sort(clusters)
+                       checkTrue( sum( srt[1:clustSize] == 1 ) / clustSize >= 0.8 )
+                       checkTrue( sum( srt[(clustSize+1):(2*clustSize)] == 2 ) / clustSize >= 0.8 )
+                       checkTrue( sum( srt[(2*clustSize+1):(3*clustSize)] == 3 ) / clustSize >= 0.8 )
+               }
+       }
+}
+
+#test several parameters agencements on custom "two moons one circle" dataset (2D)
+test.clustering3 = function()
+{
+       clustSize = 150
+
+       set.seed(32)    
+       M = matrix(nrow=3*clustSize,ncol=2)
+       #big circle: radius = 10
+       rdata = runif(clustSize, min=0, max=2*pi)
+       M[1:clustSize,] = 10 * cbind(cos(rdata), sin(rdata))
+       #moon 1: half circle of radius 5 centered at (-2, -0.5)
+       rdata = runif(clustSize, min=0, max=pi)
+       M[(clustSize+1):(2*clustSize),] = cbind(5*cos(rdata)-2, 5*sin(rdata)-0.5)
+       #moon 2: half circle of radius 5 centered at (2, 0.5)
+       rdata = runif(clustSize, min=pi, max=2*pi)
+       M[(2*clustSize+1):(3*clustSize),] = cbind(5*cos(rdata)+2, 5*sin(rdata)+0.5)
+       
+       #add small global noise
+       M = M + rnorm(2*clustSize,sd=0.1)
+       
+       #get neighborhoods from data [25 is high, but shouldn't be much lower to have 1 connex comp.]
+       NI = .Call("getNeighbors", M, 25, 0.0, 1, TRUE)
+       
+       #only ECTD distance can be used, because forcing connexity implies 
+       #creating shortcuts in graph, which strongly affect spath distance
+       distances = synclust:::getDistances("ectd", NI)
+       
+       #only hierarchical clustering can succeed here
+       clusters = synclust:::getClusters(distances, "HC", K=3)
+
+       srt = sort(clusters)
+       checkTrue( sum( srt[1:clustSize] == 1 ) / clustSize >= 0.90 )
+       checkTrue( sum( srt[(clustSize+1):(2*clustSize)] == 2 ) / clustSize >= 0.90 )
+       checkTrue( sum( srt[(2*clustSize+1):(3*clustSize)] == 3 ) / clustSize >= 0.90 )
+}
+
+#renumbering if clusters have too high labels
+test.reordering = function()
+{
+       clusters = c(1,6,8,8,8,1,1,1,6,6,6,8,8,1,1,6,8)
+       checkEquals(sort(unique(synclust:::reordering(clusters))),c(1,2,3))
+       clusters = c(23,3,23,77,77,77,1,12,12,12,77,12,23,23,12,23,77,12,23,77,1)
+       checkEquals(sort(unique(synclust:::reordering(clusters))),c(1,2,3,4,5))
+}

diff --git a/R/tests/t.connexity.R b/R/tests/t.connexity.R
new file mode 100644 (file)
index 0000000..ccabdbf
--- /dev/null
+++ b/R/tests/t.connexity.R
@@ -0,0 +1,68 @@
+#bipartite graph
+test.connexity1 = function()
+{
+       n = 10
+       NIix = as.list(rep(NA,n))
+       #connect 0 with 1, 2 with 3 ...
+       for (i in 2*(0:(n/2-1)) + 1)
+       {
+               NIix[[i]] = i+1
+               NIix[[i+1]] = i
+       }
+       cc = synclust:::getConnectedComponents(NIix)
+       #cc should contain exactly n/2 integers
+       checkEquals(n/2, length(unique(cc)))
+}
+
+#cyclic graph
+test.connexity2 = function()
+{
+       n = 10
+       NIix = as.list(rep(NA,n))
+       #connect 0 with 1, 1 with 2 ...
+       for (i in 1:n)
+               NIix[[i]] = c(ifelse(i==1,n,i-1), i%%n+1)
+       cc = synclust:::getConnectedComponents(NIix)
+       #cc should contain only one integer (1)
+       checkEquals(1, length(unique(cc)))
+}
+
+#custom graph with 3 connex components
+test.connexity3 = function()
+{
+       n = 10
+       NIix = as.list(rep(0,n))
+       NIix[[1]] = c(3,5)
+       NIix[[2]] = c(3,5)
+       NIix[[3]] = c(1,2)
+       NIix[[4]] = c(6,9,10)
+       NIix[[5]] = c(1,2)
+       NIix[[6]] = c(4)
+       NIix[[7]] = c(8)
+       NIix[[8]] = c(7)
+       NIix[[9]] = c(4)
+       NIix[[10]] = c(4,9)
+       cc = synclust:::getConnectedComponents(NIix)
+       #cc should contain only three integers
+       checkEquals(3, length(unique(cc)))
+}
+
+#custom graph, 1 connex component
+test.connexity4 = function()
+{
+       n = 10
+       NIix = as.list(rep(0,n))
+       NIix[[1]] = c(3,4,8)
+       NIix[[2]] = c(3,5,7)
+       NIix[[3]] = c(1,2)
+       NIix[[4]] = c(1,6,9,10)
+       NIix[[5]] = c(2)
+       NIix[[6]] = c(4,8)
+       NIix[[7]] = c(2)
+       NIix[[8]] = c(1,6,10)
+       NIix[[9]] = c(4)
+       NIix[[10]] = c(4,8)
+       cc = synclust:::getConnectedComponents(NIix)
+       #cc should contain only one integer (1)
+       checkEquals(1, length(unique(cc)))
+}

diff --git a/R/tests/t.utils.R b/R/tests/t.utils.R
new file mode 100644 (file)
index 0000000..3aeae6f
--- /dev/null
+++ b/R/tests/t.utils.R
@@ -0,0 +1,73 @@
+#test matrix [de]standardization
+test.de_standardize = function()
+{
+       n = 100
+       m = 10
+       M = matrix(rnorm(n*m,mean=2.0,sd=2.0),nrow=n,ncol=m)
+       
+       std = synclust:::standardize(M)
+       #result is centered:
+       checkEquals(rep(0.0, m), colMeans(std$M))
+       #result is standardized:
+       checkEquals(rep(1.0, m), sqrt( unlist( apply(std$M, 2, var) ) ))
+       
+       #rebuilt M == M:
+       M_rec = synclust:::destandardize(std)
+       checkEquals(M, M_rec)
+}
+
+#test neighborhoods remapping into one smaller component
+test.remap_neighbors = function()
+{
+       #connex comp : 1-2-5-8-10
+       #to remap into 1-2-3-4-5
+       NI = list(
+               "ix" = list(
+                       c(2,8,10), #V(1)
+                       c(1,5,8,10), #V(2)
+                       c(4,7,11), #V(3)
+                       c(3,6,9,12), #V(4)
+                       c(2,10), #V(5)
+                       c(4,7), #V(6)
+                       c(3,6,9,12), #V(7)
+                       c(1,2,10), #V(8)
+                       c(4,7,11), #V(9)
+                       c(1,2,5,8), #V(10)
+                       c(3,9), #V(11)
+                       c(4,7)), #V(12)
+               "ds" = list(
+                       c(1.0,2.0,3.0), #1
+                       c(1.0,2.0,3.0,4.0), #2
+                       c(1.0,1.0,1.0), #3
+                       c(1.0,1.0,1.0,1.0), #4
+                       c(2.0,2.0), #5
+                       c(1.0,1.0), #6
+                       c(1.0,1.0,1.0,1.0), #7
+                       c(2.0,3.0,1.0), #8
+                       c(1.0,1.0,1.0), #9
+                       c(3.0,4.0,2.0,1.0), #10
+                       c(1.0,1.0), #11
+                       c(1.0,1.0))) #12
+
+       indices = c(1,2,5,8,10)
+       locNI = synclust:::remapNeighbors(NI, indices)
+       checkEquals(2, length(locNI))
+       checkEquals(length(indices), length(locNI$ix))
+       checkEquals(length(indices), length(locNI$ds))
+       
+       #formerly index 1 (now 1)
+       checkEquals(c(2,4,5), locNI$ix[[1]])
+       checkEquals(NI$ds[[1]], locNI$ds[[1]],)
+       #formerly index 2 (now 2)
+       checkEquals(c(1,3,4,5), locNI$ix[[2]])
+       checkEquals(NI$ds[[2]], locNI$ds[[2]])
+       #formerly index 5 (now 3)
+       checkEquals(c(2,5), locNI$ix[[3]])
+       checkEquals(NI$ds[[5]], locNI$ds[[3]])
+       #formerly index 8 (now 4)
+       checkEquals(c(1,2,5), locNI$ix[[4]])
+       checkEquals(NI$ds[[8]], locNI$ds[[4]])
+       #formerly index 10 (now 5)
+       checkEquals(c(1,2,3,4), locNI$ix[[5]])
+       checkEquals(NI$ds[[10]], locNI$ds[[5]])
+}

diff --git a/R/zzz.R b/R/zzz.R
new file mode 100644 (file)
index 0000000..189dd8f
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,5 @@
+#called when package is detached ( detach("package:pkg_name") )
+.Last.lib = function(path)
+{
+       library.dynam.unload("synclust", path)
+}

diff --git a/README b/README
new file mode 100755 (executable)
index 0000000..4a85bac
--- /dev/null
+++ b/README
@@ -0,0 +1,16 @@
+Author : Benjamin Auder, Christophe Giraud
+Institution : Université Paris Sud (Orsay) / CNRS
+
+This package allow to find regions of synchronous variation for a given population species. 
+Two methods are available, either direct graph clustering or parameters estimations in a convex optimization problem.
+
+-------------------------------------------
+
+Main function to cluster populations data = findSyncVarRegions()
+
+To regenerate documentation, use R package roxygen2
+
+==================
+
+Acknowledgements :
+TODO

diff --git a/art%3A10.1007%2Fs10651-012-0222-3.pdf b/art%3A10.1007%2Fs10651-012-0222-3.pdf
new file mode 100644 (file)
index 0000000..be05319
Binary files /dev/null and b/art%3A10.1007%2Fs10651-012-0222-3.pdf differ

diff --git a/data/TODO_artificial_data b/data/TODO_artificial_data
new file mode 100644 (file)
index 0000000..e69de29

diff --git a/data/example.RData b/data/example.RData
new file mode 100644 (file)
index 0000000..7bd0146
Binary files /dev/null and b/data/example.RData differ

diff --git a/inst/TODO_manual.Rnw b/inst/TODO_manual.Rnw
new file mode 100644 (file)
index 0000000..7029b2a
--- /dev/null
+++ b/inst/TODO_manual.Rnw
@@ -0,0 +1,4 @@
+%\VignetteIndexEntry{User manual}
+\documentclass{article}
+\begin{document}
+\end{document}

diff --git a/inst/doc/convex_optimization.pdf b/inst/doc/convex_optimization.pdf
new file mode 100644 (file)
index 0000000..69e4829
Binary files /dev/null and b/inst/doc/convex_optimization.pdf differ

diff --git a/inst/doc/convex_optimization.tex b/inst/doc/convex_optimization.tex
new file mode 100644 (file)
index 0000000..7cf7a08
--- /dev/null
+++ b/inst/doc/convex_optimization.tex
@@ -0,0 +1,83 @@
+\documentclass{article}
+\usepackage{a4wide}
+\usepackage{graphicx}
+\def\pen{\textrm{pen}}
+\def\L{\mathcal{L}}
+\title{\bf Detecting areas with synchronous temporal dynamics}
+\author{Christophe Giraud}
+
+\begin{document}
+\maketitle
+
+\noindent This document summarizes the algorithm used when function \textit{findSyncVarRegions()} is called with first argument method=``convex''. Reading first the article \emph{Delimiting synchronous populations from monitoring data} by Giraud et al. is recommanded, since we use here the same notations.
+
+\section{Model and estimation procedure}
+
+\subsection{Goal}
+
+We write $Z_{stk}$ for the $k$th observations, year $t$, site $s$ and $z_{st}=\sum_{k}Z_{stk}$.
+Our goal is to estimate regions $R$ such that
+\begin{equation}\label{model}
+Z_{stk}\sim \textrm{Poisson}(\exp(\theta_{s}+f(x_{s},t)))\quad\textrm{with } f(x,t)\approx \sum_{R}\rho_{R}(t){\bf 1}_{x\in R}.
+\end{equation}
+ In other words, we try to estimate $f$ with the a priori that
+\begin{itemize}
+\item for each year $t$ the map $x \to f(x,t)$ is piecewise constant
+\item the boundary of the regions where $x \to f(x,t)$ is constant are the same for all year $t$.
+\end{itemize}
+The main difficulty is to detect the regions $R$.
+
+\subsection{Estimation procedure}
+Let $G$ be a graph and write $V(s)$ for the set of the neighbors of $s$ in G.
+The estimators $\widehat \theta$ and $\widehat f$ are defined as minimizers of 
+$$\mathcal{L}(\theta,f)+\alpha \pen(f):=\sum_{s,t}[e^{\theta_{s}+f_{st}}-z_{st}(\theta_{s}+f_{st})]+\alpha
+\sum_{s\stackrel{G}{\sim}u}\|f_{s.}-f_{u.}\|/D_{su}$$
+with boundary conditions: $f_{s1}=0$ for all $s$. We typically choose $D_{su}=1/|V(s)|+1/|V(u)|$.
+ 
+\section{Optimization algorithm}
+
+The following quantity is to be minimized
+$$\mathcal{L}(\theta,f)+\alpha \pen(f):=\sum_{s,t}[e^{\theta_{s}+f_{st}}-z_{st}(\theta_{s}+f_{st})]+\alpha\sum_{s\stackrel{G}{\sim}u}\|f_{s.}-f_{u.}\|/D_{su}$$
+with boundary conditions: $f_{s1}=0$ for all $s$.
+This last expression can be rewritten into
+$$\mathcal{L}(\theta,f)+\alpha \pen(f)=\sum_{s,t}[e^{\theta_{s}+f_{st}}-z_{st}(\theta_{s}+f_{st})]+\alpha
+\sum_{s\stackrel{G}{\sim}u}\max_{\|\phi_{su}\|\leq 1}\langle\phi_{su},f_{s.}-f_{u.}\rangle/D_{su}$$
+with $\phi_{su}\in\mathbf R^T$.
+
+\newpage
+\noindent Let us introduce
+$$F(\theta,f,\phi)=\sum_{s,t}[e^{\theta_{s}+f_{st}}-z_{st}(\theta_{s}+f_{st})]+\alpha
+\sum_{s<u}{\bf 1}_{s\stackrel{G}{\sim}u}\ \langle\phi_{su},f_{s.}-f_{u.}\rangle/D_{su}.$$
+We can reformulate the quantity to be optimized using $F$ as follows.
+$$\mathcal{L}(\theta,f)+\alpha \pen(f)=\max_{\max_{s< u}\|\phi_{su}\|\leq 1}F(\theta,f,\phi).$$
+The penalized log-likelihood can now be minimized with the following steps.
+
+\subsection*{Application}
+
+Iterate until convergence:
+\begin{enumerate}
+\item gradient descent in $\theta$:\\* $\theta\leftarrow \theta - h \nabla_{\theta}F$
+\item gradient descent in $f$ with condition $f[\ ,1]=0$\\*
+$f[\ ,-1]\leftarrow f[\ ,-1]-h'\nabla_{f[\ ,-1]}F$ 
+\item gradient ascent in $\phi$\\*
+$\phi_{su}\leftarrow \phi_{su}+h''\nabla_{\phi_{su}}F$
+\item $\phi_{su}\leftarrow \phi_{su}/\max(1,\|\phi_{su}\|)$
+\end{enumerate}
+Return($\theta,f$)
+
+%\subsection*{Gradient en $\theta$:}
+%On a
+%$$\mathcal{L}(\theta,f)=\sum_{s}\left[e^{\theta_{s}}\sum_{t}e^{f_{st}}-\theta_{s}\sum_{t}z_{st}\right]+\ldots$$
+%Donc 
+%$$\partial_{\theta_{s}}F=e^{\theta_{s}}\sum_{t}e^{f_{st}}-\sum_{t}z_{st}$$
+%
+%
+%\subsection*{Gradient en $f$:} on note $\phi_{su}=-\phi_{us}$ pour $s>u$
+%$$\partial_{f_{st}}F=e^{\theta_{s}}e^{f_{st}}-z_{st}+\alpha\sum_{u\in V(s)}\phi_{su}/D_{su}$$
+%
+%
+%\subsection*{Gradient en $\lambda$:}
+%pour $s<u$ avec $s\sim u$
+%$$\nabla_{\phi_{su}}F=\alpha(f_{s.}-f_{u.})/D_{su}.$$
+
+\end{document}

diff --git a/man/.gitignore b/man/.gitignore
new file mode 100644 (file)
index 0000000..d6b7ef3
--- /dev/null
+++ b/man/.gitignore
@@ -0,0 +1,2 @@
+*
+!.gitignore

diff --git a/src/Makevars b/src/Makevars
new file mode 100644 (file)
index 0000000..52dbd1f
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,7 @@
+PKG_CFLAGS=-g -I.
+
+PKG_LIBS=-lm -lcgds
+
+SOURCES = $(wildcard adapters/*.c sources/*.c sources/utils/*.c)
+
+OBJECTS = $(SOURCES:.c=.o)

diff --git a/src/adapters/a.connexity.c b/src/adapters/a.connexity.c
new file mode 100644 (file)
index 0000000..4c2e42e
--- /dev/null
+++ b/src/adapters/a.connexity.c
@@ -0,0 +1,49 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/connexity.h"
+
+// explore the connectivity of a graph (NIix = neighborhoods indices)
+SEXP getConnectedComponents(
+       SEXP NIix_
+) {
+       // extract NIix list vectors in a jagged array
+       int n = LENGTH(NIix_);
+       int* lengthNIix = (int*)malloc(n*sizeof(int));
+       int** NIix = (int**)malloc(n*sizeof(int*));
+       for (int i=0; i<n; i++)
+       {
+               lengthNIix[i] = LENGTH(VECTOR_ELT(NIix_,i));
+               SEXP tmp;
+               PROTECT(tmp = AS_INTEGER(VECTOR_ELT(NIix_,i)));
+               NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int));
+               for (int j=0; j<lengthNIix[i]; j++)
+                       NIix[i][j] = INTEGER(tmp)[j];
+               UNPROTECT(1);
+               // WARNING: R indices start at 1,
+               // so we must lower every index right now to avoid future bug
+               for (int j=0; j<lengthNIix[i]; j++)
+                       NIix[i][j]--;
+       }
+
+       // Main call (no R libraries)
+       int* connexComps = getConnectedComponents_core(NIix, lengthNIix, n);
+
+       // free memory
+       for (int i=0; i<n; i++)
+               free(NIix[i]);
+       free(NIix);
+       free(lengthNIix);
+
+       // transfer result in an R object
+       SEXP cc;
+       PROTECT(cc = NEW_INTEGER(n));
+       int* pcc = INTEGER_POINTER(cc);
+       for (int i=0; i<n; i++)
+               pcc[i] = connexComps[i];
+
+       // free remaining memory
+       free(connexComps);
+       UNPROTECT(1);
+
+       return cc;
+}

diff --git a/src/adapters/a.convexSolver.c b/src/adapters/a.convexSolver.c
new file mode 100644 (file)
index 0000000..9b30495
--- /dev/null
+++ b/src/adapters/a.convexSolver.c
@@ -0,0 +1,96 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/convexSolver.h"
+#include "sources/utils/algebra.h"
+
+// compute estimated ("repaired", "smoothed"...) variations from rows of M
+// NOTE: geographic coordinates dropped here, since they are unused
+SEXP getVarsWithConvexOptim(
+       SEXP M_, 
+       SEXP NIix_, 
+       SEXP alpha_, 
+       SEXP h_, 
+       SEXP epsilon_, 
+       SEXP maxiter_, 
+       SEXP symmNeighbs_, 
+       SEXP trace_
+) {
+       // get parameters
+       double alpha = NUMERIC_VALUE(alpha_);
+       double h = NUMERIC_VALUE(h_);
+       double epsilon = NUMERIC_VALUE(epsilon_);
+       int maxiter = INTEGER_VALUE(maxiter_);
+       int symmNeighbs = LOGICAL_VALUE(symmNeighbs_);
+       int trace = LOGICAL_VALUE(trace_);
+
+       // extract infos from M and get associate pointer
+       SEXP dim = getAttrib(M_, R_DimSymbol);
+       int nrow = INTEGER(dim)[0];
+       int ncol = INTEGER(dim)[1];
+       // M is always given by columns: easier to process in rows
+       double* pM = transpose(REAL(M_), nrow, ncol);
+
+       // extract NIix list vectors in a jagged array
+       int* lengthNIix = (int*)malloc(nrow*sizeof(int));
+       int** NIix = (int**)malloc(nrow*sizeof(int*));
+       for (int i=0; i<nrow; i++)
+       {
+               lengthNIix[i] = LENGTH(VECTOR_ELT(NIix_,i));
+               SEXP tmp;
+               PROTECT(tmp = AS_INTEGER(VECTOR_ELT(NIix_,i)));
+               NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int));
+               for (int j=0; j<lengthNIix[i]; j++)
+                       NIix[i][j] = INTEGER(tmp)[j];
+               UNPROTECT(1);
+               // WARNING: R indices start at 1,
+               // so we must lower every index right now to avoid future bug
+               for (int j=0; j<lengthNIix[i]; j++)
+                       NIix[i][j]--;
+       }
+
+       // Main call to core algorithm
+       Parameters params = getVarsWithConvexOptim_core(
+               pM, lengthNIix, NIix, nrow, ncol, alpha, h, epsilon, maxiter, symmNeighbs, trace);
+
+       // free neighborhoods parameters arrays
+       free(lengthNIix);
+       for (int i=0; i<nrow; i++)
+               free(NIix[i]);
+       free(NIix);
+
+       // copy matrix F into pF for output to R (1D matrices)
+       SEXP f;
+       PROTECT(f = allocMatrix(REALSXP, nrow, ncol));
+       double* pF = REAL(f);
+       for (int i=0; i<nrow; i++)
+       {
+               for (int j=0; j<ncol; j++)
+                       pF[i+nrow*j] = params.f[i][j];
+       }
+       // copy theta into pTheta for output to R
+       SEXP theta;
+       PROTECT(theta = allocVector(REALSXP, nrow));
+       double* pTheta = REAL(theta);
+       for (int i=0; i<nrow; i++)
+               pTheta[i] = params.theta[i];
+
+       // free params.f and params.theta
+       free(params.theta);
+       for (int i=0; i<nrow; i++)
+               free(params.f[i]);
+       free(params.f);
+
+       // build return list with f and theta
+       SEXP listParams, listNames;
+       PROTECT(listParams = allocVector(VECSXP, 2));
+       char* lnames[2] = {"f", "theta"}; //lists labels
+       PROTECT(listNames = allocVector(STRSXP,2));
+       for (int i=0; i<2; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+       setAttrib(listParams, R_NamesSymbol, listNames);
+       SET_VECTOR_ELT(listParams, 0, f);
+       SET_VECTOR_ELT(listParams, 1, theta);
+
+       UNPROTECT(4);
+       return listParams;
+}

diff --git a/src/adapters/a.dijkstra.c b/src/adapters/a.dijkstra.c
new file mode 100644 (file)
index 0000000..d93e591
--- /dev/null
+++ b/src/adapters/a.dijkstra.c
@@ -0,0 +1,32 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/dijkstra.h"
+
+// Dijkstra from index start : return vector of distances to every other vertex
+// NOTE [space VS perf]: passing neighborhoods infos would be enough, but
+//                                        require extra computation to translate R list argument
+SEXP dijkstra(
+       SEXP distances_, 
+       SEXP start_
+) {
+       // get arguments
+       SEXP dim = getAttrib(distances_, R_DimSymbol);
+       int n = INTEGER(dim)[0];
+       double* pDistsIn = REAL(distances_);
+       int start = INTEGER_VALUE(start_) - 1; // R indices start at 1
+
+       // Main call to core algorithm
+       double* distances = dijkstra_core(pDistsIn, start, n);
+
+       // allocate vector output and obtain R vector object
+       SEXP distsOut;
+       PROTECT(distsOut = allocVector(REALSXP, n));
+       double* pDistsOut = NUMERIC_POINTER(distsOut);
+       for (int i=0; i<n; i++)
+               pDistsOut[i] = distances[i];
+
+       free(distances);
+       UNPROTECT(1);
+
+       return distsOut;
+}

diff --git a/src/adapters/a.kmeansClustering.c b/src/adapters/a.kmeansClustering.c
new file mode 100644 (file)
index 0000000..6dc150e
--- /dev/null
+++ b/src/adapters/a.kmeansClustering.c
@@ -0,0 +1,37 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/kmeansClustering.h"
+#include <cgds/Vector.h>
+
+// k-means based on a distance matrix (nstart=10, maxiter=100)
+SEXP kmeansWithDistances(
+       SEXP distances_, 
+       SEXP K_, 
+       SEXP nstart_, 
+       SEXP maxiter_
+) {
+       // get scalar arguments
+       int K = INTEGER_VALUE(K_);
+       int nstart = NUMERIC_VALUE(nstart_);
+       int maxiter = INTEGER_VALUE(maxiter_);
+
+       // extract infos from M and get associate pointer
+       SEXP dim = getAttrib(distances_, R_DimSymbol);
+       int n = INTEGER(dim)[0];
+       double* pDistances = REAL(distances_);
+
+       // Main call to core algorithm
+       int* clusters = kmeansWithDistances_core(pDistances, n, K, nstart, maxiter);
+
+       // allocations and recopies to R vector object
+       SEXP bestClusts;
+       PROTECT(bestClusts = allocVector(INTSXP, n));
+       int* pBestClusts = INTEGER(bestClusts);
+       for (int i=0; i<n; i++)
+               pBestClusts[i] = clusters[i] + 1; // add 1 to start labels at 1
+       free(clusters);
+
+       // and return clusters
+       UNPROTECT(1);
+       return bestClusts;
+}

diff --git a/src/adapters/a.neighbors.c b/src/adapters/a.neighbors.c
new file mode 100644 (file)
index 0000000..1de503c
--- /dev/null
+++ b/src/adapters/a.neighbors.c
@@ -0,0 +1,79 @@
+#include <R.h>
+#include <Rdefines.h>
+#include "sources/neighbors.h"
+#include "sources/utils/algebra.h"
+#include <cgds/List.h>
+
+// Function to obtain neighborhoods.
+// NOTE: alpha = weight parameter to compute distances; -1 means "adaptive"
+// WARNING : M is given in columns
+SEXP getNeighbors(
+       SEXP M_, 
+       SEXP k_, 
+       SEXP alpha_, 
+       SEXP gmode_, 
+       SEXP simpleDists_
+) {
+       // get scalar arguments
+       int k = INTEGER_VALUE(k_);
+       double alpha = NUMERIC_VALUE(alpha_);
+       int gmode = INTEGER_VALUE(gmode_);
+       int simpleDists = LOGICAL_VALUE(simpleDists_);
+
+       // extract infos from M and get associate pointer
+       SEXP dim = getAttrib(M_, R_DimSymbol);
+       int nrow = INTEGER(dim)[0];
+       int ncol = INTEGER(dim)[1];
+       // M is always given by columns: easier to process in rows
+       double* pM = transpose(REAL(M_), nrow, ncol);
+
+       // Main call to core algorithm which fills neighborhoods lists
+       List** neighborhoods = getNeighbors_core(pM, alpha, k, gmode, simpleDists, nrow, ncol);
+
+       // transfer neighborhoods lists into R vectors
+       SEXP NIix, NIds;
+       PROTECT(NIix = allocVector(VECSXP, nrow)); //indices
+       PROTECT(NIds = allocVector(VECSXP, nrow)); //distances
+       for (int i=0; i<nrow; i++)
+       {
+               SEXP neighbsIX, neighbsDS;
+               PROTECT(neighbsIX = NEW_INTEGER(list_size(neighborhoods[i])));
+               PROTECT(neighbsDS = NEW_NUMERIC(list_size(neighborhoods[i])));
+               int* pNeighbsIX = INTEGER_POINTER(neighbsIX);
+               double* pNeighbsDS = NUMERIC_POINTER(neighbsDS);
+               ListIterator* neighbsI = list_get_iterator(neighborhoods[i]);
+               int j = 0;
+               while (listI_has_data(neighbsI))
+               {
+                       IndDist indDist; listI_get(neighbsI, indDist);
+                       // WARNING: R arrays start at index 1
+                       pNeighbsIX[j] = indDist.index + 1;
+                       pNeighbsDS[j] = indDist.distance;
+                       j++;
+                       listI_move_next(neighbsI);
+               }
+               SET_VECTOR_ELT(NIix, i, neighbsIX);
+               SET_VECTOR_ELT(NIds, i, neighbsDS);
+               UNPROTECT(2);
+               listI_destroy(neighbsI);
+               list_destroy(neighborhoods[i]);
+       }
+       free(neighborhoods);
+
+       // create R list labels to access with NI$ix and NI$ds
+       SEXP listNames;
+       char* lnames[2] = {"ix", "ds"}; //lists labels
+       PROTECT(listNames = allocVector(STRSXP,2));
+       for (int i=0; i<2; i++)
+               SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
+
+       // allocate and fill neighborhoods list to return
+       SEXP NI;
+       PROTECT(NI = allocVector(VECSXP, 2));
+       SET_VECTOR_ELT(NI, 0, NIix);
+       SET_VECTOR_ELT(NI, 1, NIds);
+       setAttrib(NI, R_NamesSymbol, listNames);
+
+       UNPROTECT(4);
+       return NI;
+}

diff --git a/src/sources/connexity.c b/src/sources/connexity.c
new file mode 100644 (file)
index 0000000..db9a185
--- /dev/null
+++ b/src/sources/connexity.c
@@ -0,0 +1,53 @@
+#include "sources/connexity.h"
+#include <cgds/Stack.h>
+#include <stdlib.h>
+#include "sources/utils/boolean.h"
+
+int* getConnectedComponents_core(int** NIix, int* lengthNIix, int n)
+{
+       int* cc = (int*)calloc(n,sizeof(int));
+       Stack* toBeExplored = stack_new(int);
+       int curInd = 0, nextInd;
+       bool* alreadyExpanded = (bool*)calloc(n,sizeof(bool));
+
+       // while the entire graph hasn't been visited
+       while (S_TRUE)
+       {
+               int label = curInd+1;
+               cc[curInd] = label;
+               stack_push(toBeExplored, curInd);
+
+               // while new elements are discovered in current component,
+               // mark them as expanded and stack their neighbors
+               while (!stack_empty(toBeExplored))
+               {
+                       stack_top(toBeExplored, nextInd);
+                       stack_pop(toBeExplored);
+                       cc[nextInd] = label;
+
+                       for (int j=0; j<lengthNIix[nextInd]; j++)
+                       {
+                               if (!alreadyExpanded[NIix[nextInd][j]])
+                                       stack_push(toBeExplored, NIix[nextInd][j]);
+                       }
+                       alreadyExpanded[nextInd] = S_TRUE;
+               }
+
+               // curInd is set to the next unexplored index (if possible)
+               for (int i=0; i<n; i++)
+               {
+                       if (cc[i] == 0)
+                       {
+                               curInd = i;
+                               break;
+                       }
+               }
+               if (cc[curInd] != 0)
+                       break;
+       }
+
+       free(alreadyExpanded);
+       stack_destroy(toBeExplored);
+
+       return cc;
+}

diff --git a/src/sources/connexity.h b/src/sources/connexity.h
new file mode 100644 (file)
index 0000000..6660d73
--- /dev/null
+++ b/src/sources/connexity.h
@@ -0,0 +1,11 @@
+#ifndef SYNCLUST_CONNEXITY_H
+#define SYNCLUST_CONNEXITY_H
+
+// Core algorithm to find connected components in an undirected graph
+int* getConnectedComponents_core(
+       int** NIix, 
+       int* lengthNIix, 
+       int n
+);
+
+#endif

diff --git a/src/sources/convexSolver.c b/src/sources/convexSolver.c
new file mode 100644 (file)
index 0000000..0661e08
--- /dev/null
+++ b/src/sources/convexSolver.c
@@ -0,0 +1,219 @@
+#include "sources/convexSolver.h"
+#include <stdio.h> //to trace LL evolution
+#include <stdlib.h>
+#include <math.h>
+#include "sources/utils/algebra.h"
+
+// auxiliary to compute log-likelihood + penalty
+double computeLogLikelihood(
+       double** f, double* theta, double** Z, double*** phi, 
+       int* lengthNIix, int** NIix, double alpha, int nrow, int ncol) 
+{
+       double LL = 0.0;
+
+       // for each row in data matrix: (one row = observations from 2001 to 2009)
+       for (int i=0; i<nrow; i++)
+       {
+               // theta[i] == -INFINITY if no birds were seen at this site
+               if (theta[i] != -INFINITY)
+               {
+                       // for each year from 2001 to 2009:
+                       for (int j=0; j<ncol; j++)
+                               LL += (exp(theta[i] + f[i][j]) - Z[i][j] * (theta[i] + f[i][j]));
+               }
+               // add penalty term
+               double penalty = 0.0;
+               double Ds = 1.0/lengthNIix[i];
+               // lengthNIix[i] == size of the neighborhood of site i
+               for (int j=0; j<lengthNIix[i]; j++)
+               {
+                       // compute <phi[s,u] , f[s,] - f[u,]> with u == NIix[i][j] (j-th neighbor of i)
+                       double dotProduct = 0.0;
+                       for (int jj=0; jj<ncol; jj++)
+                               dotProduct += phi[i][NIix[i][j]][jj] * (f[i][jj] - f[NIix[i][j]][jj]);
+                       // normalization by sum of inverses of neighborhoods sizes
+                       double Dsu = Ds + 1.0/lengthNIix[NIix[i][j]];
+                       penalty += dotProduct / Dsu;
+               }
+               LL += alpha * penalty;
+       }
+
+       return LL;
+}
+
+// compute estimated ("repaired", "smoothed"...) variations from rows of M
+// NOTE: geographic coordinates dropped here, since they are unused
+Parameters getVarsWithConvexOptim_core(
+                       double* pM, int* lengthNIix, int** NIix, int nrow, int ncol, 
+                       double alpha, double h, double epsilon, int maxiter, bool symmNeighbs, bool trace)
+{
+       double EPS = 1e-10; // HACK: some small numerical constant to avoid oddities
+
+       // theta_s = log(average z_st)
+       double* theta = (double*)calloc(nrow,sizeof(double));
+       // NOTE:: Z is 'double' because it is [can be] an average value (of integers)
+       double** Z = (double**)malloc(nrow*sizeof(double*));
+       for (int i=0; i<nrow; i++)
+       {
+               Z[i] = (double*)malloc(ncol*sizeof(double));
+               for (int j=0; j<ncol; j++)
+               {
+                       Z[i][j] = pM[i*ncol+j];
+                       theta[i] += Z[i][j];
+               }
+               // since pM values are assumed to be integers (and ncol not too high ?!),
+               // the following test may be simplified into (theta[i]==0.0)
+               if (fabs(theta[i]) < EPS)
+                       theta[i] = -INFINITY;
+               else
+                       theta[i] = log(theta[i]/ncol);
+       }
+       // initialize f to observed variations
+       double** F = (double**)malloc(nrow*sizeof(double*));
+       for (int i=0; i<nrow; i++)
+       {
+               F[i] = (double*)calloc(ncol,sizeof(double));
+               if (theta[i] != -INFINITY)
+               {
+                       for (int j=0; j<ncol; j++)
+                       {
+                               if (Z[i][j] > 0.0)
+                                       F[i][j] = log(Z[i][j]) - theta[i];
+                       }
+               }
+       }
+       // phi_s,u = 1/sqrt(ncol) (1 ... 1) [TODO:: costly in memory !]
+       double invSqrtNcol = 1.0/sqrt(ncol);
+       double*** phi = (double***)malloc(nrow*sizeof(double**));
+       for (int i=0; i<nrow; i++)
+       {
+               phi[i] = (double**)malloc(nrow*sizeof(double*));
+               for (int ii=0; ii<nrow; ii++)
+               {
+                       phi[i][ii] = (double*)malloc(ncol*sizeof(double));
+                       for (int j=0; j<ncol; j++)
+                               phi[i][ii][j] = invSqrtNcol;
+               }
+       }
+
+       // initialize log-likelihood
+       double LL = computeLogLikelihood(
+               F, theta, Z, phi, lengthNIix, NIix, alpha, nrow, ncol);
+       double oldLL = -INFINITY;
+
+       /*******************
+        * START ITERATIONS
+        *******************/
+
+       int kounter = 0; // limit iterations count, in case of
+       while (fabs(LL - oldLL) >= epsilon && kounter++ < maxiter)
+       {
+               // gradient descent for theta
+               for (int i=0; i<nrow; i++) {
+                       if (theta[i] == -INFINITY)
+                               continue; // skip these sites: cannot get information
+                       double sumExpFst = 0.0;
+                       for (int j=0; j<ncol; j++)
+                               sumExpFst += exp(F[i][j]);
+                       double sumZst = 0.0;
+                       for (int j=0; j<ncol; j++)
+                               sumZst += Z[i][j];
+                       double gradI = exp(theta[i]) * sumExpFst - sumZst;
+                       theta[i] -= h * gradI;
+               }
+
+               // gradient descent for f
+               double sumDdivPhi;
+               for (int i=0; i<nrow; i++)
+               {
+                       double invDs = 1.0/lengthNIix[i];
+                       for (int j=0; j<ncol; j++)
+                       {
+                               double gradIJ = - Z[i][j];
+                               if (theta[i] != -INFINITY)
+                               {
+                                       // no theta[i] contribution if nothing observed
+                                       gradIJ += exp(theta[i] + F[i][j]);
+                               }
+                               // + sum on neighbors u: s-->u, - sum on neighbors u: u-->s
+                               sumDdivPhi = 0.0;
+                               for (int jj=0; jj<lengthNIix[i]; jj++)
+                               {
+                                       double Dsu = invDs + 1.0/lengthNIix[NIix[i][jj]];
+                                       sumDdivPhi += phi[i][NIix[i][jj]][j] / Dsu;
+                                       if (symmNeighbs)
+                                       {
+                                               //shortcut: if symmetric neighborhoods, it's easy to sum on u-->s
+                                               sumDdivPhi -= phi[NIix[i][jj]][i][j] / Dsu;
+                                       }
+                               }
+                               gradIJ += alpha * sumDdivPhi;
+                               if (!symmNeighbs)
+                               {
+                                       // need to remove sum on neighbors u: u-->s, uneasy way.
+                                       //TODO: computation is much too costly here; need preprocessing
+                                       sumDdivPhi = 0.0;
+                                       for (int ii=0; ii<nrow; ii++)
+                                       {
+                                               //~ if (ii == i) continue;
+                                               for (int jj=0; jj<lengthNIix[ii]; jj++)
+                                               {
+                                                       if (NIix[ii][jj] == i)
+                                                       {
+                                                               sumDdivPhi += phi[ii][i][j] / (invDs + 1.0/lengthNIix[ii]);
+                                                               break; //i can only appear once among neighbors of ii
+                                                       }
+                                               }
+                                       }
+                                       gradIJ -= alpha * sumDdivPhi;
+                               }
+                               F[i][j] -= h * gradIJ;
+                       }
+               }
+
+               // gradient ascent for phi
+               for (int i=0; i<nrow; i++)
+               {
+                       double Ds = 1.0/lengthNIix[i];
+                       for (int ii=0; ii<nrow; ii++)
+                       {
+                               double Dsu = Ds + 1.0/lengthNIix[ii];
+                               for (int j=0; j<ncol; j++)
+                               {
+                                       double gradI_II_J = alpha * (F[i][j] - F[ii][j]) / Dsu;
+                                       phi[i][ii][j] += h * gradI_II_J;
+                               }
+                               // also renormalize to have ||phi_su|| == 1.0
+                               double normPhi = norm2(phi[i][ii], ncol);
+                               //~ if (normPhi > 1.0) {
+                               if (normPhi > EPS)
+                               {
+                                       for (int j=0; j<ncol; j++)
+                                               phi[i][ii][j] /= normPhi;
+                               }
+                       }
+               }
+
+               oldLL = LL;
+               LL = computeLogLikelihood(
+                       F, theta, Z, phi, lengthNIix, NIix, alpha, nrow, ncol);
+               if (trace)
+                       printf("%i / LLs: %.0f %.0f\n",kounter,oldLL,LL); // optional trace of LL evolution
+       } /*** END ITERATIONS ***/
+
+       // free all local parameters arrays but (theta, F) (used as return value)
+       for (int i=0; i<nrow; i++)
+       {
+               free(Z[i]);
+               for (int ii=0; ii<nrow; ii++)
+                       free(phi[i][ii]);
+               free(phi[i]);
+       }
+       free(Z);
+       free(phi);
+
+       Parameters params;
+       params.f = F;
+       params.theta = theta;
+       return params;
+}

diff --git a/src/sources/convexSolver.h b/src/sources/convexSolver.h
new file mode 100644 (file)
index 0000000..ac59a3a
--- /dev/null
+++ b/src/sources/convexSolver.h
@@ -0,0 +1,54 @@
+#ifndef SYNCLUST_CONVEXSOLVER_H
+#define SYNCLUST_CONVEXSOLVER_H
+
+#include "sources/utils/boolean.h"
+
+// auxiliary to compute euclidian norm
+double norm2(
+       double* v, 
+       int length
+);
+
+// auxiliary to compute euclidian distance
+double distance2(
+       double* f1, 
+       double* f2, 
+       int length
+);
+
+// auxiliary to compute log-likelihood + penalty
+double computeLogLikelihood(
+       double** f, 
+       double* theta, 
+       double** Zst, 
+       double*** phi, 
+       int* lengthNIix, 
+       int** NIix, 
+       double alpha, 
+       int nrow, 
+       int ncol
+);
+
+// structure to return parameters (theta, f) [and others if needed later]
+typedef struct Parameters {
+    double** f;
+    double* theta;
+} Parameters;
+
+// compute estimated ("repaired", "smoothed"...) variations from rows of M
+// NOTE: geographic coordinates dropped here, since they are unused
+Parameters getVarsWithConvexOptim_core(
+       double* pM, 
+       int* lengthNIix, 
+       int** NIix, 
+       int nrow, 
+       int ncol, 
+       double alpha, 
+       double h, 
+       double epsilon, 
+       int maxiter, 
+       bool symmNeighbs, 
+       bool trace
+);
+
+#endif

diff --git a/src/sources/dijkstra.c b/src/sources/dijkstra.c
new file mode 100644 (file)
index 0000000..4207c36
--- /dev/null
+++ b/src/sources/dijkstra.c
@@ -0,0 +1,55 @@
+#include "sources/dijkstra.h"
+#include <stdlib.h>
+#include "sources/utils/boolean.h"
+#include <math.h>
+
+// Dijkstra from index start : return vector of distances to every other vertex
+// TODO: use a good priority queue, and pass NI instead of pDistsIn (+ linear preprocessing)
+double* dijkstra_core(double* pDistsIn, int start, int n) {
+
+       // initalisations
+       double* pDistsOut = (double*)malloc(n*sizeof(double));
+       bool* visited = (bool*)malloc(n*sizeof(bool));
+       for (int i=0; i<n; i++) {
+               pDistsOut[i] = INFINITY;
+               visited[i] = S_FALSE; // nothing seen so far
+       }
+       pDistsOut[start] = 0.0; // we are at distance 0 from self
+
+       while (S_TRUE)
+       {
+               double minGeodDist = INFINITY;
+
+               // n1 <-- node in "unseen" with smallest geodesic distance
+               // NOTE: on first loop, n1 == start
+               int n1 = 0;
+               for (int i=0; i<n; i++)
+               {
+                       if (!visited[i] && pDistsOut[i] < minGeodDist)
+                       {
+                               n1 = i;
+                               minGeodDist = pDistsOut[i];
+                       }
+               }
+
+               if (minGeodDist == INFINITY)
+                       break; // all is explored
+
+               visited[n1] = S_TRUE; // n1 is expanded
+
+               // For n2 running through neighbors of n1
+               for (int n2 = 0; n2<n; n2++)
+               {
+                       int ind_n12 = n1*n+n2; // or n1+n*n2 (symmetry)
+                       if (!isnan(pDistsIn[ind_n12]))
+                       {
+                               // check if we'd better go through n1 (on the way from start to n2)
+                               if (pDistsOut[n2] > pDistsOut[n1] + pDistsIn[ind_n12])
+                                       pDistsOut[n2] = pDistsOut[n1] + pDistsIn[ind_n12];
+                       }
+               }
+       }
+       free(visited);
+
+       return pDistsOut;
+}

diff --git a/src/sources/dijkstra.h b/src/sources/dijkstra.h
new file mode 100644 (file)
index 0000000..fbcedaa
--- /dev/null
+++ b/src/sources/dijkstra.h
@@ -0,0 +1,11 @@
+#ifndef SYNCLUST_DIJKSTRA_H
+#define SYNCLUST_DIJKSTRA_H
+
+// Dijkstra from index start : return vector of distances to every other vertex
+double* dijkstra_core(
+       double* pDistsIn, 
+       int start, 
+       int n
+);
+
+#endif

diff --git a/src/sources/kmeansClustering.c b/src/sources/kmeansClustering.c
new file mode 100644 (file)
index 0000000..c7dec16
--- /dev/null
+++ b/src/sources/kmeansClustering.c
@@ -0,0 +1,199 @@
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+#include "sources/utils/boolean.h"
+#include "sources/kmeansClustering.h"
+
+// auxiliary function to obtain a random sample of 1..n with K elements
+void sample(int* centers, int n, int K)
+{
+       // refVect = (0,1,...,n-1,n)
+       int* refVect = (int*)malloc(n*sizeof(int));
+       for (int i=0; i<n; i++)
+               refVect[i] = i;
+
+       int curSize = n; // current size of the sampling set
+       for (int j=0; j<K; j++)
+       {
+               // pick an index in sampling set:
+               int index = rand()%curSize;
+               centers[j] = refVect[index];
+               // move this index outside of sampling set:
+               refVect[index] = refVect[--curSize];
+       }
+
+       free(refVect);
+}
+
+// auxiliary function to compare two sets of centers
+int unequalCenters(int* ctrs1, int* ctrs2, int n, int K)
+{
+       // HACK: special case at initialization, ctrs2 = 0
+       if (K > 1 && ctrs2[0]==0 && ctrs2[1]==0)
+               return S_TRUE;
+
+       // compVect[i] == 1 iff index i is found in ctrs1 or ctrs2
+       int* compVect = (int*)calloc(n,sizeof(int));
+
+       int kountNonZero = 0; // count non-zero elements in compVect
+       for (int j=0; j<K; j++)
+       {
+               if (compVect[ctrs1[j]] == 0)
+                       kountNonZero++;
+               compVect[ctrs1[j]] = 1;
+               if (compVect[ctrs2[j]] == 0)
+                       kountNonZero++;
+               compVect[ctrs2[j]] = 1;
+       }
+
+       free(compVect);
+
+       // if we found more than K non-zero elements, ctrs1 and ctrs2 differ
+       return (kountNonZero > K);
+}
+
+// assign a vector (represented by its distances to others, as distances[index,])
+// to a cluster, represented by its center ==> output is integer in 0..K-1
+int assignCluster(int index, double* distances, int* centers, int n, int K)
+{
+       int minIndex = 0;
+       double minDist = distances[index*n+centers[0]];
+
+       for (int j=1; j<K; j++)
+       {
+               if (distances[index*n+centers[j]] < minDist)
+               {
+                       minDist = distances[index*n+centers[j]];
+                       minIndex = j;
+               }
+       }
+
+       return minIndex;
+}
+
+// k-means based on a distance matrix (nstart=10, maxiter=100)
+int* kmeansWithDistances_core(
+       double* distances, int n, int K, int nstart, int maxiter)
+{
+       int* bestClusts = (int*)malloc(n*sizeof(int));
+       double bestDistor = INFINITY;
+       double avgClustSize = (double)n/K;
+       int* ctrs = (int*)malloc(K*sizeof(int));
+       int* oldCtrs = (int*)malloc(K*sizeof(int));
+       Vector** clusters = (Vector**)malloc(K*sizeof(Vector*));
+       for (int j=0; j<K; j++)
+               clusters[j] = vector_new(int);
+
+       // set random number generator seed
+       srand(time(NULL));
+
+       for (int startKount=0; startKount < nstart; startKount++)
+       {
+               // centers (random) [re]initialization
+               sample(ctrs,n,K);
+               for (int j=0; j<K; j++)
+                       oldCtrs[j] = 0;
+               int kounter = 0;
+
+               /*************
+                *  main loop
+                *************/
+
+               // while old and "new" centers differ..
+               while (unequalCenters(ctrs,oldCtrs,n,K) && kounter++ < maxiter)
+               {
+                       // (re)initialize clusters to empty sets
+                       for (int j=0; j<K; j++)
+                               vector_clear(clusters[j]);
+
+                       // estimate clusters belongings
+                       for (int i=0; i<n; i++)
+                       {
+                               int affectation = assignCluster(i, distances, ctrs, n, K);
+                               vector_push(clusters[affectation], i);
+                       }
+
+                       // copy current centers to old centers
+                       for (int j=0; j<K; j++)
+                               oldCtrs[j] = ctrs[j];
+
+                       // recompute centers
+                       for (int j=0; j<K; j++)
+                       {
+                               int minIndex = -1;
+                               double minSumDist = INFINITY;
+                               VectorIterator* iter1 = vector_get_iterator(clusters[j]);
+                               vectorI_reset_begin(iter1);
+                               while (vectorI_has_data(iter1))
+                               {
+                                       int index1; vectorI_get(iter1, index1);
+                                       // attempt to use current index as center
+                                       double sumDist = 0.0;
+                                       VectorIterator* iter2 = vector_get_iterator(clusters[j]);
+                                       vectorI_reset_begin(iter2);
+                                       while (vectorI_has_data(iter2))
+                                       {
+                                               int index2; vectorI_get(iter2, index2);
+                                               sumDist += distances[index1*n+index2];
+                                               vectorI_move_next(iter2);
+                                       }
+                                       if (sumDist < minSumDist)
+                                       {
+                                               minSumDist = sumDist;
+                                               minIndex = index1;
+                                       }
+                                       vectorI_destroy(iter2);
+                                       vectorI_move_next(iter1);
+                               }
+                               if (minIndex >= 0)
+                                       ctrs[j] = minIndex;
+                               // HACK: some 'random' index (a cluster should never be empty)
+                               // this case should never happen anyway 
+                               // (pathological dataset with replicates)
+                               else
+                                       ctrs[j] = 0;
+                               vectorI_destroy(iter1);
+                       }
+               } /***** end main loop *****/
+
+               // finally compute distorsions :
+               double distor = 0.0;
+               for (int j=0; j<K; j++)
+               {
+                       VectorIterator* iter = vector_get_iterator(clusters[j]);
+                       vectorI_reset_begin(iter);
+                       while (vectorI_has_data(iter))
+                       {
+                               int index; vectorI_get(iter, index);
+                               distor += distances[index*n+ctrs[j]];
+                               vectorI_move_next(iter);
+                       }
+                       vectorI_destroy(iter);
+               }
+               if (distor < bestDistor)
+               {
+                       // copy current clusters into bestClusts
+                       for (int j=0; j<K; j++)
+                       {
+                               VectorIterator* iter = vector_get_iterator(clusters[j]);
+                               vectorI_reset_begin(iter);
+                               while (vectorI_has_data(iter))
+                               {
+                                       int index; vectorI_get(iter, index);
+                                       bestClusts[index] = j;
+                                       vectorI_move_next(iter);
+                               }
+                               vectorI_destroy(iter);
+                       }
+                       bestDistor = distor;
+               }
+       }
+
+       free(ctrs);
+       free(oldCtrs);
+       for (int j=0; j<K; j++)
+               vector_destroy(clusters[j]);
+       free(clusters);
+
+       return bestClusts;
+}

diff --git a/src/sources/kmeansClustering.h b/src/sources/kmeansClustering.h
new file mode 100644 (file)
index 0000000..e2fd346
--- /dev/null
+++ b/src/sources/kmeansClustering.h
@@ -0,0 +1,40 @@
+#ifndef SYNCLUST_KMEANSCLUSTERING_H
+#define SYNCLUST_KMEANSCLUSTERING_H
+
+#include <cgds/Vector.h>
+
+// auxiliary function to obtain a random sample of 1..n with K elements
+void sample(
+       int* centers, 
+       int n, 
+       int K
+);
+
+// auxiliary function to compare two sets of centers
+int unequalCenters(
+       int* ctrs1, 
+       int* ctrs2, 
+       int n, 
+       int K
+);
+
+// assign a vector (represented by its distances to others, as distances[index,])
+// to a cluster, represented by its center index ==> output is integer in 0..K-1
+int assignCluster(
+       int index, 
+       double* distances, 
+       int* centers, 
+       int n, 
+       int K
+);
+
+// k-means based on a distance matrix (nstart=10, maxiter=100)
+int* kmeansWithDistances_core(
+       double* pDistances, 
+       int n, 
+       int K, 
+       int nstart, 
+       int maxiter
+);
+
+#endif

diff --git a/src/sources/neighbors.c b/src/sources/neighbors.c
new file mode 100644 (file)
index 0000000..d9407a4
--- /dev/null
+++ b/src/sources/neighbors.c
@@ -0,0 +1,216 @@
+#include "sources/neighbors.h"
+#include <stdlib.h>
+#include <math.h>
+#include "sources/utils/algebra.h"
+#include <cgds/BufferTop.h>
+#include "sources/utils/boolean.h"
+
+// evaluate distance between M[i,] and M[ii,]
+double getDistance(double* M, int i, int ii, int ncol, double alpha, 
+       bool simpleDists)
+{
+       // if simpleDists is ON, it means we are in stage 2 after convex optimization
+       // ==> use full data, we know that now there are no NA's.
+       if (simpleDists)
+               return distance2(M+i*ncol, M+ii*ncol, ncol);
+
+       // get distance for values per year
+       double dist1 = 0.0;
+       int valCount = 0; // number of not-NA fields
+       int nobs = ncol-2; // ncol is 9+2 for our initial dataset (2001 to 2009)
+       for (int j=0; j<nobs; j++)
+       {
+               if (!isnan(M[i*ncol+j]) && !isnan(M[ii*ncol+j]))
+               {
+                       double diff = M[i*ncol+j] - M[ii*ncol+j];
+                       dist1 += diff*diff;
+                       valCount++;
+               }
+       }
+       if (valCount > 0)
+               dist1 /= valCount;
+
+       // get distance for coordinates values
+       double dist2 = 0.0;
+       for (int j=nobs; j<ncol; j++)
+       {
+               double diff = M[i*ncol+j] - M[ii*ncol+j];
+               dist2 += diff*diff;
+       }
+       dist2 /= 2.0; //to harmonize with above normalization
+       if (valCount == 0)
+               return sqrt(dist2); // no other choice
+
+       //NOTE: adaptive alpha, the more NA's in vector, 
+       //        the more geo. coords. are taken into account
+       alpha = (alpha >= 0.0 ? alpha : (double)valCount/nobs);
+       return sqrt(alpha*dist1 + (1.0-alpha)*dist2);
+}
+
+// symmetrize neighborhoods lists (augmenting or reducing)
+void symmetrizeNeighbors(List** neighborhoods, int nrow, int gmode)
+{
+       IndDist curNeighbI, curNeighbJ;
+       for (int i=0; i<nrow; i++)
+       {
+               ListIterator* iterI = list_get_iterator(neighborhoods[i]);
+               while (listI_has_data(iterI))
+               {
+                       listI_get(iterI, curNeighbI);
+                       // check if neighborhoods[curNeighbI->index] has i
+                       bool reciproc = S_FALSE;
+                       List* neighbsJ = neighborhoods[curNeighbI.index];
+                       ListIterator* iterJ = list_get_iterator(neighbsJ);
+                       while (listI_has_data(iterJ))
+                       {
+                               listI_get(iterJ, curNeighbJ);
+                               if (curNeighbJ.index == i)
+                               {
+                                       reciproc = S_TRUE;
+                                       break;
+                               }
+                               listI_move_next(iterJ);
+                       }
+
+                       if (!reciproc)
+                       {
+                               if (gmode == 1)
+                               {
+                                       // augmenting:
+                                       // add (previously) non-mutual neighbor to neighbsJ
+                                       list_insert_back(neighbsJ, i);
+                               }
+                               // test list_size() >= 2 because we don't allow empty neighborhoods
+                               else if (gmode == 0 && list_size(neighborhoods[i]) >= 2)
+                               {
+                                       // reducing:
+                                       // remove non-mutual neighbor to neighbsI
+                                       listI_remove(iterI,BACKWARD);
+                               }
+                       }
+                       listI_move_next(iterI);
+                       listI_destroy(iterJ);
+               }
+               listI_destroy(iterI);
+       }
+}
+
+// restrain neighborhoods: choose one per quadrant (for convex optimization)
+void restrainToQuadrants(List** neighborhoods, int nrow, int ncol, double* M)
+{
+       IndDist curNeighbI;
+       for (int i=0; i<nrow; i++)
+       {
+               ListIterator* iter = list_get_iterator(neighborhoods[i]);
+               // choose one neighbor in each quadrant (if available); 
+               // WARNING: multi-constraint optimization,
+               //   > as close as possible to angle bissectrice
+               //   > not too far from current data point
+
+               // resp. SW,NW,SE,NE "best" neighbors :
+               int bestIndexInDir[4] = {-1,-1,-1,-1};
+               // corresponding "performances" :
+               double bestPerfInDir[4] = {INFINITY,INFINITY,INFINITY,INFINITY};
+               while (listI_has_data(iter))
+               {
+                       listI_get(iter, curNeighbI);
+                       // get delta_x and delta_y to know in which quadrant 
+                       // we are and then get "index performance"
+                       // ASSUMPTION: all sites are distinct
+                       double deltaX = 
+                               M[i*ncol+(ncol-2)] - M[curNeighbI.index*ncol+(ncol-2)];
+                       double deltaY = 
+                               M[i*ncol+(ncol-1)] - M[curNeighbI.index*ncol+(ncol-1)];
+                       double angle = fabs(atan(deltaY/deltaX));
+                       // naive combination; [TODO: improve]
+                       double perf = curNeighbI.distance + fabs(angle-M_PI_4);
+                       // map {-1,-1} to 0, {-1,1} to 1 ...etc :
+                       int index = 2*(deltaX>0)+(deltaY>0);
+                       if (perf < bestPerfInDir[index])
+                       {
+                               bestIndexInDir[index] = curNeighbI.index;
+                               bestPerfInDir[index] = perf;
+                       }
+                       listI_move_next(iter);
+               }
+
+               // restrain neighborhood to the "best directions" found
+               listI_reset_head(iter);
+               while (listI_has_data(iter))
+               {
+                       listI_get(iter, curNeighbI);
+                       // test list_size() <= 1 because we don't allow empty neighborhoods
+                       if (list_size(neighborhoods[i]) <= 1 ||
+                               curNeighbI.index==bestIndexInDir[0] || 
+                               curNeighbI.index==bestIndexInDir[1] || 
+                               curNeighbI.index==bestIndexInDir[2] || 
+                               curNeighbI.index==bestIndexInDir[3])
+                       {
+                               // OK, keep it
+                               listI_move_next(iter);
+                               continue;
+                       }
+                       // remove current node
+                       listI_remove(iter,FORWARD);
+               }
+               listI_destroy(iter);
+       }
+}
+
+// Function to obtain neighborhoods.
+// NOTE: alpha = weight parameter to compute distances; -1 means "adaptive"
+List** getNeighbors_core(double* M, double alpha, int k, int gmode, 
+       bool simpleDists, int nrow, int ncol)
+{
+       // prepare list buffers to get neighborhoods
+       // (OK for small to moderate values of k)
+       BufferTop** bufferNeighbs = 
+               (BufferTop**)malloc(nrow*sizeof(BufferTop*));
+       for (int i=0; i<nrow; i++)
+               bufferNeighbs[i] = buffertop_new(IndDist, k, MIN_T, 2);
+
+       // MAIN LOOP
+
+       // for each row in M, find its k nearest neighbors
+       for (int i=0; i<nrow; i++)
+       {
+               // for each potential neighbor...
+               for (int ii=0; ii<nrow; ii++)
+               {
+                       if (ii == i) 
+                               continue;
+
+                       // evaluate distance from M[i,] to M[ii,]
+                       double distance = 
+                               getDistance(M, i, ii, ncol, alpha, simpleDists);
+
+                       // (try to) add index 'ii' + distance to bufferNeighbs[i]
+                       IndDist id = (IndDist){.index=ii, .distance=distance};
+                       buffertop_tryadd(bufferNeighbs[i], id, distance);
+               }
+       }
+
+       // free buffers and transfer their contents into lists easier to process
+       List** neighborhoods = (List**)malloc(nrow*sizeof(List*));
+       for (int i=0; i<nrow; i++)
+       {
+               neighborhoods[i] = buffertop_2list(bufferNeighbs[i]);
+               buffertop_destroy(bufferNeighbs[i]);
+       }
+       free(bufferNeighbs);
+
+       // OPTIONAL MUTUAL KNN
+       if (gmode==0 || gmode==1)
+       {
+               // additional processing to symmetrize neighborhoods (augment or not)
+               symmetrizeNeighbors(neighborhoods, nrow, gmode);
+       }
+       else if (gmode==3)
+       {
+               // choose one neighbor per quadrant (for convex optimization)
+               restrainToQuadrants(neighborhoods, nrow, ncol, M);
+       }
+       // nothing to do if gmode==2 (simple assymmetric kNN)
+
+       return neighborhoods;
+}

diff --git a/src/sources/neighbors.h b/src/sources/neighbors.h
new file mode 100644 (file)
index 0000000..9cd39f3
--- /dev/null
+++ b/src/sources/neighbors.h
@@ -0,0 +1,50 @@
+#ifndef SYNCLUST_NEIGHBORS_H
+#define SYNCLUST_NEIGHBORS_H
+
+#include "sources/utils/boolean.h"
+#include <cgds/List.h>
+
+// evaluate distance between M[i,] and M[ii,]
+double getDistance(
+       double* M, 
+       int i, 
+       int ii, 
+       int ncol, 
+       double alpha, 
+       bool simpleDists
+);
+
+// symmetrize neighborhoods lists (augmenting or reducing)
+void symmetrizeNeighbors(
+       List** neighborhoods, 
+       int nrow, 
+       int gmode
+);
+
+// restrain neighborhoods: choose one per quadrant (for convex optimization)
+void restrainToQuadrants(
+       List** neighborhoods, 
+       int nrow, 
+       int ncol, 
+       double* M
+);
+
+// structure to store a neighbor index and the distance to this neighbor
+typedef struct IndDist {
+       int index;
+       double distance;
+} IndDist;
+
+// Function to obtain neighborhoods.
+// NOTE: alpha = weight parameter to compute distances; -1 means "adaptive"
+List** getNeighbors_core(
+       double* M, 
+       double alpha, 
+       int k, 
+       int gmode, 
+       bool simpleDists, 
+       int nrow, 
+       int ncol
+);
+
+#endif

diff --git a/src/sources/utils/algebra.c b/src/sources/utils/algebra.c
new file mode 100644 (file)
index 0000000..4a28273
--- /dev/null
+++ b/src/sources/utils/algebra.c
@@ -0,0 +1,37 @@
+#include "sources/utils/algebra.h"
+#include <stdlib.h>
+#include <math.h>
+
+// small useful function to transform a matrix as given by R 
+// into a easier-to-handle one.
+double* transpose(double* M, int nrow, int ncol)
+{
+       double* Mtr = (double*)malloc(nrow*ncol*sizeof(double));
+       for (int i=0; i<nrow; i++)
+       {
+               for (int j=0; j<ncol; j++)
+                       Mtr[i*ncol+j] = M[i+nrow*j];
+       }
+       return Mtr;
+}
+
+// auxiliary to compute euclidian norm
+double norm2(double* v, int length)
+{
+       double norm = 0.0;
+       for (int j=0; j<length; j++)
+               norm += v[j]*v[j];
+       return sqrt(norm);
+}
+
+// auxiliary to compute euclidian distance
+double distance2(double* v1, double* v2, int length)
+{
+       double distance = 0.0, diff;
+       for (int j=0; j<length; j++)
+       {
+               diff = v1[j]-v2[j];
+               distance += diff*diff;
+       }
+       return sqrt(distance);
+}

diff --git a/src/sources/utils/algebra.h b/src/sources/utils/algebra.h
new file mode 100644 (file)
index 0000000..c64a9dc
--- /dev/null
+++ b/src/sources/utils/algebra.h
@@ -0,0 +1,25 @@
+#ifndef SYNCLUST_ALGEBRA_H
+#define SYNCLUST_ALGEBRA_H
+
+// small useful function to transform a matrix as given by R 
+// into a easier-to-handle one.
+double* transpose(
+       double* M, 
+       int nrow, 
+       int ncol
+);
+
+// auxiliary to compute euclidian norm
+double norm2(
+       double* v, 
+       int length
+);
+
+// auxiliary to compute euclidian distance
+double distance2(
+       double* v1, 
+       double* v2, 
+       int length
+);
+
+#endif

diff --git a/src/sources/utils/boolean.h b/src/sources/utils/boolean.h
new file mode 100644 (file)
index 0000000..63ac5d2
--- /dev/null
+++ b/src/sources/utils/boolean.h
@@ -0,0 +1,10 @@
+#ifndef SYNCLUST_BOOLEAN_H
+#define SYNCLUST_BOOLEAN_H
+
+// boolean type (prefix because of conflict with R booleans)
+typedef enum {
+       S_FALSE=0,
+       S_TRUE=1
+} bool;
+
+#endif

diff --git a/src/tests/.gitignore b/src/tests/.gitignore
new file mode 100644 (file)
index 0000000..23af2c7
--- /dev/null
+++ b/src/tests/.gitignore
@@ -0,0 +1 @@
+testexec

diff --git a/src/tests/Makefile b/src/tests/Makefile
new file mode 100644 (file)
index 0000000..c0d6b14
--- /dev/null
+++ b/src/tests/Makefile
@@ -0,0 +1,22 @@
+CC = gcc
+CFLAGS = -std=gnu99 -g
+LDFLAGS = -L.. -lcgds -lm ../synclust.so
+INCLUDES = -I..
+TARGET = testexec
+
+all: testexec
+
+clean:
+       rm -f *.o
+       rm -f testexec
+
+.PHONY: all clean
+
+SOURCES = main.c helpers.c t.connexity.c t.convexSolver.c t.dijkstra.c t.kmeansClustering.c t.neighbors.c
+OBJECTS = $(SOURCES:%.c=%.o)
+
+testexec: $(OBJECTS)
+       $(CC) $(LDFLAGS) -o $@ $^ 
+
+%.o: %.c
+       $(CC) $(CFLAGS) $(INCLUDES) -o $@ -c $<

diff --git a/src/tests/helpers.c b/src/tests/helpers.c
new file mode 100644 (file)
index 0000000..310c9ca
--- /dev/null
+++ b/src/tests/helpers.c
@@ -0,0 +1,110 @@
+#include "tests/helpers.h"
+
+//auxiliary to check vectors equality
+int checkEqualV(double* v1, double* v2, int n)
+{
+       double epsilon = 1e-10; // arbitrary small value to make comparisons
+       for (int i=0; i<n; i++)
+       {
+               if (fabs(v1[i] - v2[i]) >= epsilon)
+                       return S_FALSE;
+       }
+       return S_TRUE;
+}
+
+// auxiliary to count distinct values in an integer array
+int countDistinctValues(int* v, int n)
+{
+       int maxVal = v[0];
+       for (int i=1; i<n; i++)
+       {
+               if (v[i] > maxVal)
+                       maxVal = v[i];
+       }
+       int* kountArray = (int*)calloc(maxVal+1,sizeof(int));
+       int res = 0;
+       for (int i=0; i<n; i++)
+       {
+               if (kountArray[v[i]] == 0)
+               {
+                       res++;
+                       kountArray[v[i]] = 1; // mark this value as "seen"
+               }
+       }
+       free(kountArray);
+       return res;
+}
+
+// helper to check clustering result
+int labelIsProcessed(int label, int* processedLabels, int countProcessedLabels)
+{
+       for (int j=0; j<countProcessedLabels; j++)
+       {
+               if (processedLabels[j] == label)
+                       return S_TRUE;
+       }
+       return S_FALSE;
+}
+
+// check both clusters purity and size (ideally they all contain n/clustCount points)
+int checkClustersProportions(int* clusters, int n, int clustCount, double tol)
+{
+       // initialize array to keep track of clusters labels
+       int* processedLabels = (int*)malloc(clustCount*sizeof(int));
+       for (int j=0; j<clustCount; j++)
+               processedLabels[j] = -1;
+       int countProcessedLabels = 0, clustSize = n/clustCount;
+
+       int i=0;
+       while (i<n)
+       {
+               // go to the next unprocessed label (if possible)
+               while (i < n && labelIsProcessed(clusters[i], processedLabels, countProcessedLabels))
+               {
+                       i++;
+               }
+               if (i >= n)
+                       break;
+
+               int label = clusters[i];
+               processedLabels[countProcessedLabels++] = label;
+
+               // count elements in current cluster (represented by label)
+               int countDataWithCurLabel = 0;
+               for (int ii=0; ii<n; ii++)
+               {
+                       if (clusters[ii] == label)
+                               countDataWithCurLabel++;
+               }
+               // if cluster is too small or too big, test fails
+               if ( fabs(countDataWithCurLabel - clustSize) / n > tol )
+               {
+                       free(processedLabels);
+                       return S_FALSE;
+               }
+
+               // now check counts per cluster (repartition);
+               // the labels should not be spread between different (true) groups
+               int maxCountLabelInClust = 0;
+               for (int kounter=0; kounter<clustCount; kounter++)
+               {
+                       int countLabelInClust = 0;
+                       for (int ii=kounter*clustSize; ii<(kounter+1)*clustSize; ii++)
+                       {
+                               if (clusters[ii] == label)
+                                       countLabelInClust++;
+                       }
+                       if (countLabelInClust > maxCountLabelInClust)
+                               maxCountLabelInClust = countLabelInClust;
+               }
+               // we must have max(repartition) / clustSize >= 1 - tol
+               if ((double)maxCountLabelInClust / clustSize < 1.0 - tol)
+               {
+                       free(processedLabels);
+                       return S_FALSE;
+               }
+       }
+
+       free(processedLabels);
+       return S_TRUE;
+}

diff --git a/src/tests/helpers.h b/src/tests/helpers.h
new file mode 100644 (file)
index 0000000..3d3f5aa
--- /dev/null
+++ b/src/tests/helpers.h
@@ -0,0 +1,49 @@
+#ifndef SYNCLUST_HELPERS_H
+#define SYNCLUST_HELPERS_H
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#include "sources/utils/boolean.h"
+
+//auxiliary to check vectors equality
+int checkEqualV(
+       double* v1, 
+       double* v2, 
+       int n
+);
+
+// auxiliary to count distinct values in an integer array
+int countDistinctValues(
+       int* v, 
+       int n
+);
+
+// check if clusters proportions match a given fraction (tolerance 'tol')
+int checkClustersProportions(
+       int* clusters, 
+       int n, 
+       int clustCount, 
+       double tol
+);
+
+// basic unit tests macros
+#define ASSERT_TRUE(b) \
+do { \
+       if (!(b)) { \
+               fprintf(stdout, "Error in file %s at line %i\n", __FILE__, __LINE__); \
+               return; \
+       } \
+} while (0)
+
+#define ASSERT_FALSE(b) \
+do { \
+       if (b) { \
+               fprintf(stdout, "Error in file %s at line %i\n", __FILE__, __LINE__); \
+               return; \
+       } \
+} while (0)
+
+#endif

diff --git a/src/tests/main.c b/src/tests/main.c
new file mode 100644 (file)
index 0000000..beb9ed6
--- /dev/null
+++ b/src/tests/main.c
@@ -0,0 +1,27 @@
+void connexityTests();
+void convexSolverTests();
+void dijkstraTests();
+void kmeansClusteringTests();
+void neighborsTests();
+
+// Main function calling all unit tests
+// If nothing gets printed, everything's OK
+int main()
+{
+       test_connexity1();
+       test_connexity2();
+
+       test_convexSolver1();
+       test_convexSolver2();
+
+       test_dijkstra1();
+       test_dijkstra2();
+
+       test_kmeansClustering1();
+       test_kmeansClustering2();
+
+       test_neighbors1();
+       test_neighbors2();
+
+       return 0;
+}

diff --git a/src/tests/t.connexity.c b/src/tests/t.connexity.c
new file mode 100644 (file)
index 0000000..ae5e98f
--- /dev/null
+++ b/src/tests/t.connexity.c
@@ -0,0 +1,117 @@
+#include "tests/helpers.h"
+#include "sources/connexity.h"
+
+//completely disconnected graph (no edges)
+void test_connexity1()
+{
+       int n = 10;
+       int* lengthNIix = (int*)calloc(n,sizeof(int));
+       int** NIix = (int**)malloc(n*sizeof(int*));
+       for (int i=0; i<n; i++) NIix[i] = NULL;
+       int* cc = getConnectedComponents_core(NIix, lengthNIix, n);
+       //cc should contain all integers from 1 to n
+       ASSERT_TRUE(countDistinctValues(cc,n) == n);
+       free(lengthNIix);
+       free(NIix);
+       free(cc);
+}
+
+//bipartite graph
+void test_connexity2()
+{
+       int n = 10;
+       int* lengthNIix = (int*)malloc(n*sizeof(int));
+       int** NIix = (int**)malloc(n*sizeof(int*));
+       for (int i=0; i<n; i++)
+       {
+               lengthNIix[i] = 1;
+               NIix[i] = (int*)malloc(sizeof(int));
+       }
+       // connect 0 with 1, 2 with 3 ...
+       for (int i=0; i<n; i+=2)
+       {
+               NIix[i][0] = i+1;
+               NIix[i+1][0] = i;
+       }
+       int* cc = getConnectedComponents_core(NIix, lengthNIix, n);
+       //cc should contain exactly n/2 integers
+       ASSERT_TRUE(countDistinctValues(cc,n) == n/2);
+       free(lengthNIix);
+       for (int i=0; i<n; i++)
+               free(NIix[i]);
+       free(NIix);
+       free(cc);
+}
+
+//~ //cyclic graph
+//~ void test_connexity3() {
+       //~ int n = 10;
+       //~ int* adjMat = (int*)calloc(n*n,sizeof(int));
+       //~ // connect 0 with 1, 1 with 2 ...
+       //~ for (int i=0; i<n; i++) {
+               //~ adjMat[i+n*((i+1)%n)] = TRUE;
+               //~ adjMat[(i+1)%n+n*i] = TRUE;
+       //~ }
+       //~ int* cc = getConnectedComponents_core(adjMat, n);
+       //~ //cc should contain only one value (presumably '1')
+       //~ warnIfFails(countDistinctValues(cc,n) == 1, "c3", "");
+       //~ free(adjMat);
+       //~ free(cc);
+//~ }
+//~ 
+//~ //custom graph with 3 connex components
+//~ void test_connexity4() {
+       //~ int n = 10;
+       //~ 
+       //~ 
+               //~ {2,4},
+               //~ {2,4},
+               //~ {0,1},
+               //~ {5,8,9},
+               //~ {0,1},
+               //~ {3,7},
+               //~ {},
+               //~ {5,9},
+               //~ {3},
+               //~ {3,7}
+       //~ 
+       //~ int adjMat[100] = 
+       //~ {
+               //~ 0,0,1,0,1,0,0,0,0,0,
+               //~ 0,0,1,0,1,0,0,0,0,0,
+               //~ 1,1,0,0,0,0,0,0,0,0,
+               //~ 0,0,0,0,0,1,0,0,1,1,
+               //~ 1,1,0,0,0,0,0,0,0,0,
+               //~ 0,0,0,1,0,0,0,1,0,0,
+               //~ 0,0,0,0,0,0,0,0,0,0,
+               //~ 0,0,0,0,0,1,0,0,0,1,
+               //~ 0,0,0,1,0,0,0,0,0,0,
+               //~ 0,0,0,1,0,0,0,1,0,0
+       //~ };
+       //~ int* cc = getConnectedComponents_core(adjMat, n);
+       //~ //cc should contain exactly 3 values
+       //~ warnIfFails(countDistinctValues(cc,n) == 3, "c4", "");
+       //~ free(cc);
+//~ }
+//~ 
+//~ //custom graph, 1 connex component
+//~ void test_connexity5() {
+       //~ int n = 10;
+       //~ int adjMat[100] = 
+       //~ {
+               //~ 0,0,1,1,0,0,0,1,0,0,
+               //~ 0,0,1,0,1,0,1,0,0,0,
+               //~ 1,1,0,0,0,0,0,0,0,0,
+               //~ 1,0,0,0,0,1,0,0,1,1,
+               //~ 0,1,0,0,0,0,0,0,0,0,
+               //~ 0,0,0,1,0,0,0,1,0,0,
+               //~ 0,1,0,0,0,0,0,0,0,0,
+               //~ 1,0,0,0,0,1,0,0,0,1,
+               //~ 0,0,0,1,0,0,0,0,0,0,
+               //~ 0,0,0,1,0,0,0,1,0,0
+       //~ };
+       //~ int* cc = getConnectedComponents_core(adjMat, n);
+       //~ //cc should contain only one value (presumably '1')
+       //~ warnIfFails(countDistinctValues(cc,n) == 1, "c5", "");
+       //~ free(cc);
+//~ }

diff --git a/src/tests/t.convexSolver.c b/src/tests/t.convexSolver.c
new file mode 100644 (file)
index 0000000..679d78b
--- /dev/null
+++ b/src/tests/t.convexSolver.c
@@ -0,0 +1,749 @@
+#include "tests/helpers.h"
+#include "sources/convexSolver.h"
+#include "sources/neighbors.h"
+
+//test on random matrix
+void test_convexSolver1()
+{
+       int nrow=100, ncol=11;
+       //usual parameters
+       double alpha = 0.2;
+       int k = 4, gmode = 3;
+
+       // initialize M to random matrix which somewhat simulate true data repartition
+       srand(time(NULL));
+       double* M = (double*)malloc(nrow*ncol*sizeof(double));
+       // fill series first
+       for (int j=0; j<ncol-2; j++)
+       {
+               for (int i=0; i<nrow; i++)
+                       M[i+nrow*j] = (double)rand()/RAND_MAX;
+       }
+       // fill coordinates
+       int minCoord[2] = {  77100.0, 1023100.0}; // X, Y
+       int maxCoord[2] = {1710200.0, 2672200.0}; // X, Y
+       for (int j=ncol-2; j<ncol; j++)
+       {
+               for (int i=0; i<nrow; i++)
+               {
+                       M[i+nrow*j] = minCoord[ncol-j-1] + 
+                               (maxCoord[ncol-j-1]-minCoord[ncol-j-1]) * (double)rand()/RAND_MAX;
+               }
+       }
+
+       // we need neighborhoods to apply solver
+       List** NI_tmp = getNeighbors_core(M, alpha, k, gmode, S_FALSE, nrow, ncol);
+       int** NIix = (int**)malloc(nrow*sizeof(int*));
+       int* lengthNIix = (int*)malloc(nrow*sizeof(int));
+       for (int i=0; i<nrow; i++)
+       {
+               ListIterator* iterJ = list_get_iterator(NI_tmp[i]);
+               int neighbSize = list_size(NI_tmp[i]);
+               lengthNIix[i] = neighbSize;
+               NIix[i] = (int*)malloc(neighbSize*sizeof(int));
+               for (int j=0; j<neighbSize; j++)
+               {
+                       int tmpNeighb; listI_get(iterJ, tmpNeighb);
+                       NIix[i][j] = tmpNeighb;
+                       listI_move_next(iterJ);
+               }
+               listI_destroy(iterJ);
+               list_destroy(NI_tmp[i]);
+       }
+       free(NI_tmp);
+
+       // restrict M by removing coordinates columns
+       double* M_restrict = (double*)malloc(nrow*(ncol-2)*sizeof(double));
+       for (int j=0; j<ncol-2; j++)
+       {
+               for (int i=0; i<nrow; i++)
+                       M_restrict[i+nrow*j] = M[i+nrow*j];
+       }
+       free(M);
+
+       // finally apply test [TODO?]
+       double pcoef = 1.0, h = 1e-3, eps = 1e-4;
+       int maxit = 1e3;
+       Parameters params = getVarsWithConvexOptim_core(
+               M_restrict, lengthNIix, NIix, nrow, ncol-2, pcoef, h, eps, maxit, S_FALSE, S_TRUE);
+
+       // free parameters memory
+       free(M_restrict);
+       free(lengthNIix);
+       for (int i=0; i<nrow; i++)
+               free(NIix[i]);
+       free(NIix);
+
+       // free results memory
+       free(params.theta);
+       for (int i=0; i<nrow; i++)
+               free(params.f[i]);
+       free(params.f);
+}
+
+//test with real dataset (birds observations)
+void test_convexSolver2()
+{
+       int nrow=625, ncol=9;
+       //usual parameters
+       double alpha = 0.2;
+       int k = 10, gmode = 1;
+
+       // initialize M to completed data from birds observations [NO coordinates]
+       double M[] = {0.20,1.0,1.10,0.70,0.60,0.0,0.0,0.10,0.50,0.50,0.0,1.0,1.30,1.0,0.0,1.50,0.90,0.40,0.70,1.10,1.20,0.90,1.10,0.80,1.20,1.70,1.60,1.30,1.30,0.30,1.20,0.30,0.20,0.30,0.80,0.0,0.30,0.30,1.20,0.80,0.0,0.0,0.90,0.90,1.50,0.50,0.30,1.90,1.20,1.80,1.20,3.40,2.20,1.10,0.90,0.80,0.50,1.70,1.50,0.60,0.90,0.90,0.60,0.0,1.50,1.0,0.20,0.30,1.10,1.10,0.50,2.80,0.50,0.80,0.30,0.70,2.0,1.0,1.30,0.40,0.0,0.80,1.0,1.40,0.50,1.0,1.10,1.60,1.10,0.80,1.0,0.0,2.50,1.60,1.20,0.80,1.0,1.20,0.40,0.60,0.0,0.20,0.40,1.40,1.90,1.10,1.30,1.70,0.70,0.10,1.90,0.0,0.0,0.0,0.0,0.80,0.30,0.60,0.70,0.20,0.40,2.20,2.50,0.50,1.0,1.80,1.20,1.30,0.666667,0.90,1.30,0.70,1.50,0.70,0.40,1.10,0.80,0.30,1.0,1.10,2.0,1.30,1.80,0.0,1.60,0.60,0.80,1.10,0.80,2.0,1.40,2.333333,1.70,0.60,0.70,1.40,0.70,0.50,1.120,0.40,1.28750,0.08750,0.08750,0.60,1.780,1.70,2.60,1.840,0.716667,0.133333,0.30,0.0,2.0,0.144444,0.20,1.70,0.58750,0.871429,0.416667,0.30,0.614286,0.714286,0.10,1.057143,1.233333,0.90,1.433333,0.60,1.16250,2.528571,1.5250,1.36250,1.150,0.942857,1.28750,1.30,2.23750,1.70,1.316667,1.50,1.70,1.350,0.520,0.60,1.385714,1.550,0.483333,1.20,1.36250,2.30,2.0,2.7750,2.933333,2.533333,1.80,0.333333,2.21250,1.2750,1.350,1.366667,1.20,2.10,0.083333,2.3750,1.414286,1.01250,1.650,1.01250,1.23750,1.46250,0.80,0.550,0.60,0.783333,0.3250,1.10,0.620,0.80,0.483333,1.057143,0.20,0.6250,0.3250,0.033333,0.250,0.150,0.20,0.666667,1.0,0.050,0.60,1.950,0.357143,0.0750,0.871429,1.0250,1.90,0.138889,0.90,0.76250,0.90,0.350,0.840,0.183333,1.60,0.90,0.46250,0.30,0.20,1.985714,2.26250,1.040,0.7750,1.550,1.433333,1.1750,0.30,2.20,0.4750,1.20,0.950,0.850,0.314286,1.30,0.866667,0.850,0.70,1.050,1.142857,0.0,0.8750,1.40,2.20,2.250,2.0,0.10,0.0,0.828571,0.26250,0.250,1.171429,0.6250,0.4250,1.10,0.8750,1.30,1.30,0.333333,0.642857,1.1750,1.250,0.242857,1.133333,0.60,0.20,1.442857,0.31250,0.450,1.23750,0.73750,1.033333,2.20,2.3750,0.30,2.26250,1.91250,1.78750,0.50,1.16250,1.7750,2.0,1.20,1.1250,0.70,1.10,0.871429,0.950,2.814286,0.96250,0.6750,1.2750,0.5250,1.46250,1.050,0.43750,0.350,0.30,1.380,0.66250,0.483333,1.56250,1.414286,1.457143,1.33750,1.120,1.350,0.81250,0.816667,2.73750,0.80,2.91250,0.30,0.608333,1.191837,2.2250,1.98750,2.4750,2.183333,2.250,3.5750,0.614286,0.03750,0.10,0.250,0.0,1.450,1.70,1.4250,2.214286,4.60,2.5250,1.10,1.30,1.80,1.350,1.33750,2.142857,1.50,1.9750,1.750,1.085714,1.7750,1.250,0.260,1.63750,1.250,1.70,0.742857,2.11250,1.050,0.60,0.533333,0.6250,0.820,1.060,0.880,0.10,1.585714,1.250,0.633333,0.90,0.48750,0.0,0.8750,1.528571,2.750,2.060,0.68750,1.70,1.76250,2.28750,2.86250,2.23750,1.60,1.4250,0.750,0.86250,1.13750,1.48750,0.90,0.914286,0.666667,0.70,0.0,1.30,0.60,0.842857,1.742857,2.080,0.166667,0.814286,2.050,1.60,2.580,2.10,2.233333,0.985714,1.785714,1.066667,1.90,2.180,0.70,1.60,0.780,2.10,0.071429,0.70,1.066667,0.840,0.816667,0.366667,0.80,0.40,1.50,0.883333,1.433333,0.071429,0.340,0.014286,1.10,1.228571,0.90,1.233333,1.466667,0.747619,0.80,0.928571,1.30,2.014286,0.50,1.166667,0.850,0.650,0.185714,0.550,0.10,1.485714,0.485714,0.0,1.10,0.850,0.857143,0.30,0.457143,0.20,1.446032,0.50,0.80,1.10,0.966667,0.550,1.150,0.950,0.714286,1.916667,0.533333,1.866667,0.850,0.750,0.20,0.980,0.90,0.880,0.266667,1.160,0.666667,1.50,1.740,1.680,1.3250,0.666667,1.240,0.560,0.533333,0.550,0.040,1.150,0.2250,0.90,1.050,0.460,1.466667,0.10,0.350,2.483333,1.666667,0.166667,0.0,1.30,1.340,2.0,1.480,1.460,1.633333,0.70,1.020,1.1750,0.30,1.680,1.7750,1.250,0.920,1.060,0.380,1.520,0.860,0.10,1.580,0.50,1.0750,0.333333,0.40,0.0750,1.950,0.9250,0.750,0.20,1.466667,1.150,0.4750,2.60,1.633333,1.30,0.566667,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.20,0.70,0.70,0.70,0.60,0.10,0.0,0.20,0.50,0.80,0.10,0.70,0.50,0.70,0.0,0.90,0.80,0.60,0.70,0.90,1.30,1.10,1.0,0.90,1.40,2.70,1.70,1.80,2.80,0.80,0.70,0.70,0.444444,0.40,0.70,0.20,0.20,0.40,0.90,0.80,0.0,0.0,1.10,1.20,1.40,0.60,0.40,0.90,1.20,2.10,1.50,3.40,2.20,1.20,1.0,1.40,0.50,2.0,1.30,0.80,0.80,0.90,0.80,0.60,1.50,1.20,0.20,0.30,1.0,0.80,0.40,2.0,1.0,1.10,0.80,1.0,2.20,1.10,1.40,0.40,0.542857,0.80,1.10,1.40,0.50,1.0,1.10,1.10,1.20,0.50,1.228571,0.0,2.30,2.0,1.10,0.90,0.90,1.20,1.10,0.40,0.10,0.10,2.40,1.70,2.250,1.20,2.70,1.90,0.20,0.40,1.50,0.30,0.0,0.0,0.0,0.750,0.30,0.7750,0.70,0.80,0.40,2.20,2.50,0.50,1.40,1.40,1.40,1.30,0.40,0.385714,2.30,1.20,2.0,1.60,0.10,1.0,1.0,0.30,1.20,2.80,2.0,0.80,2.0,0.0,2.10,1.10,1.0,1.0,1.0,1.80,2.60,1.961905,1.80,0.70,0.850,1.0,0.70,1.40,1.20,0.20,1.20,0.0,0.0,0.60,2.70,2.20,2.40,2.30,1.0,0.0,0.30,0.0,2.30,0.10,0.20,1.0,0.20,1.0,0.30,0.20,0.30,0.60,0.10,1.10,1.50,0.80,0.30,0.20,0.60,2.90,1.40,1.80,1.20,1.40,1.10,1.0,2.80,1.90,1.10,1.50,1.60,1.30,0.80,0.30,1.10,2.30,1.10,1.20,1.80,2.10,2.40,1.80,3.20,2.90,1.90,0.30,2.30,1.10,1.30,1.50,0.90,2.40,0.0,1.80,2.10,1.30,1.80,1.0,1.30,1.50,1.10,0.30,0.70,1.10,0.10,1.0,0.70,0.80,0.70,1.30,0.20,0.70,0.0,0.0,0.10,0.10,0.10,0.40,0.90,0.0,0.90,2.50,0.10,0.0,0.80,0.50,1.60,0.10,0.90,0.70,1.10,0.30,0.80,0.0,1.70,0.90,0.60,0.20,0.20,1.40,2.60,0.30,0.60,1.60,1.40,1.30,0.40,1.40,0.60,1.90,1.60,0.60,0.30,1.30,0.90,0.90,0.60,1.0,1.0,0.0,0.90,1.40,2.20,2.60,2.10,0.0,0.0,0.50,0.50,0.20,1.0,0.80,0.60,0.90,0.90,1.30,1.40,0.70,0.90,1.30,1.0,0.30,1.50,0.50,0.20,1.90,0.10,0.50,0.60,1.0,1.10,2.20,2.10,0.20,1.60,1.50,1.60,0.50,1.40,1.70,2.10,1.20,1.0,0.70,1.30,0.60,0.80,3.20,1.30,1.0,1.30,0.60,1.50,0.80,0.0,0.30,0.30,1.0,0.60,0.50,1.70,0.60,1.40,1.10,1.10,1.10,0.90,1.10,2.30,0.90,2.30,0.30,0.70,1.10,1.90,1.70,2.20,2.60,2.30,2.80,0.40,0.0,0.30,0.30,0.0,1.90,1.80,1.30,1.60,2.30,0.80,1.10,1.20,1.30,1.30,1.50,2.10,1.60,2.10,1.60,0.50,2.0,1.0,0.40,0.80,1.10,2.30,0.70,1.30,1.10,0.20,0.40,0.40,0.70,1.0,0.70,0.0,1.90,1.20,0.70,0.90,0.60,0.0,0.60,1.70,2.60,2.0,0.30,1.50,1.30,1.80,1.80,1.90,1.60,1.10,1.10,0.90,0.80,1.0,0.90,0.50,0.666667,0.70,0.0,1.30,0.60,0.842857,1.742857,2.080,0.166667,0.814286,2.050,1.60,2.580,2.10,2.233333,0.985714,1.785714,1.066667,1.90,2.180,0.70,1.60,0.780,2.10,0.071429,0.70,1.066667,0.840,0.816667,0.366667,0.80,0.40,1.50,0.883333,1.433333,0.071429,0.340,0.014286,1.10,1.228571,0.90,1.233333,1.466667,0.747619,0.80,0.928571,1.30,2.014286,0.50,1.166667,0.850,0.650,0.185714,0.550,0.10,1.485714,0.485714,0.0,1.10,0.850,0.857143,0.30,0.457143,0.20,1.446032,0.50,0.80,1.10,0.966667,0.550,1.150,0.950,0.714286,1.916667,0.533333,1.866667,0.850,0.750,0.20,0.980,0.90,0.880,0.266667,1.160,0.666667,1.50,1.740,1.680,1.3250,0.666667,1.240,0.560,0.533333,0.550,0.040,1.150,0.2250,0.90,1.050,0.460,1.466667,0.10,0.350,2.483333,1.666667,0.166667,0.0,1.30,1.340,2.0,1.480,1.460,1.633333,0.70,1.020,1.1750,0.30,1.680,1.7750,1.250,0.920,1.060,0.380,1.520,0.860,0.10,1.580,0.50,1.0750,0.333333,0.40,0.0750,1.950,0.9250,0.750,0.20,1.466667,1.150,0.4750,2.60,1.633333,1.30,0.566667,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.10,0.90,0.50,1.0,0.80,0.050,0.10,0.150,0.50,0.650,0.10,0.850,0.50,0.70,0.0,0.80,1.40,0.50,0.50,0.80,1.70,1.40,1.80,1.60,1.20,2.10,2.20,2.80,1.70,0.80,0.90,0.20,0.40,0.50,0.80,0.20,0.10,0.60,0.90,1.10,0.0,0.0,1.0,1.050,1.30,1.40,0.60,1.0,1.40,1.950,1.50,3.40,2.20,1.0,1.0,1.98750,0.557143,2.30,1.30,0.90,0.80,1.40,0.40,0.80,1.70,1.0,0.20,0.60,0.90,1.40,0.70,2.0,0.40,1.0,0.30,1.40,1.60,1.70,1.70,0.40,0.20,0.80,0.90,1.40,0.50,1.0,1.10,1.40,1.0,0.70,0.90,0.0,2.30,2.20,0.90,1.0,0.70,1.20,0.30,0.40,0.10,0.10,2.40,2.10,2.50,1.40,2.0,1.80,0.20,0.30,1.50,0.30,0.0,0.0,0.0,0.70,0.30,1.10,0.70,0.80,0.40,2.20,2.50,0.60,1.20,1.20,0.60,1.30,0.30,0.385714,1.20,1.20,1.50,1.0,0.50,0.90,1.20,0.30,1.30,1.20,2.40,1.20,1.90,0.0,1.90,1.10,0.90,1.11250,1.0,1.60,2.0,2.0,1.60,0.70,0.850,1.10,0.70,0.40,1.0,0.20,1.40,0.10,0.10,0.60,1.60,1.70,2.70,1.50,1.20,0.10,0.10,0.0,2.30,0.0,0.20,1.90,0.80,0.70,0.40,0.20,0.40,0.50,0.20,0.80,1.30,0.90,1.60,0.70,1.40,2.50,1.90,1.50,1.10,1.10,1.70,1.80,2.0,1.70,1.30,1.50,1.70,1.40,0.40,0.50,1.40,0.80,0.70,1.20,1.60,2.30,2.50,3.40,2.60,2.10,1.90,0.60,2.40,1.60,1.40,0.90,1.30,1.90,0.0,1.90,1.30,0.90,2.0,1.0,1.10,0.80,0.90,0.80,0.60,0.90,0.50,1.10,0.50,0.80,0.70,1.70,0.20,0.70,0.10,0.10,0.20,0.40,0.30,0.80,1.10,0.10,0.30,2.30,0.60,0.20,0.70,1.10,2.30,0.0,0.60,1.0,0.90,0.40,0.60,0.50,2.20,0.90,0.60,0.20,0.20,1.90,2.60,0.40,0.90,1.550,2.10,1.50,0.60,2.20,0.30,1.50,1.10,1.10,0.10,1.30,0.80,0.80,1.20,0.40,0.90,0.0,0.60,1.40,2.20,2.30,2.50,0.0,0.0,1.0,0.20,0.10,1.30,0.60,0.30,1.40,0.60,1.30,1.40,0.10,0.80,1.0,1.30,0.30,1.10,0.60,0.20,1.50,0.50,0.40,0.80,0.60,1.20,2.20,2.20,0.30,2.0,1.60,1.70,0.50,1.10,2.60,2.50,1.20,1.1250,0.70,1.10,0.60,0.80,3.0,0.80,0.40,1.10,0.50,1.20,1.20,0.30,0.20,0.30,1.40,0.40,0.50,1.50,0.30,1.40,1.50,1.40,0.80,0.90,0.90,2.20,1.10,2.70,0.30,0.60,1.20,2.80,2.0,3.0,2.60,2.20,3.20,1.30,0.0,0.10,0.20,0.0,1.0,1.60,1.40,2.0,3.0,3.40,1.10,1.30,2.30,1.60,1.30,1.70,1.60,4.20,1.90,0.80,1.80,1.20,0.20,1.0,1.30,2.0,0.40,1.80,0.80,0.50,0.80,0.40,1.30,0.90,0.90,0.30,2.90,1.30,0.50,0.90,0.40,0.0,0.80,1.80,2.90,2.20,0.70,1.40,1.70,2.80,2.70,2.0,1.60,1.10,0.80,0.50,1.0,1.50,0.90,1.10,0.80,0.70,0.0,1.30,0.60,0.70,1.20,1.90,0.20,0.50,1.30,0.90,2.70,1.70,1.90,1.10,0.90,0.70,1.0,1.70,0.90,1.70,0.80,2.10,0.10,0.70,1.60,1.30,0.80,0.60,0.50,0.40,1.20,0.90,1.70,0.0,0.20,0.0,0.90,1.20,0.70,1.30,1.30,1.10,0.40,0.70,1.0,1.70,0.40,0.80,0.80,0.80,0.20,0.20,0.0,1.50,0.70,0.0,1.0,1.20,0.40,0.30,0.0,0.20,1.10,0.50,0.50,1.30,0.90,0.40,1.50,0.90,1.30,2.20,0.10,1.866667,0.850,0.750,0.20,0.980,0.90,0.880,0.266667,1.160,0.666667,1.50,1.740,1.680,1.3250,0.666667,1.240,0.560,0.533333,0.550,0.040,1.150,0.2250,0.90,1.050,0.460,1.466667,0.10,0.350,2.483333,1.666667,0.166667,0.0,1.30,1.340,2.0,1.480,1.460,1.633333,0.70,1.020,1.1750,0.30,1.680,1.7750,1.250,0.920,1.060,0.380,1.520,0.860,0.10,1.580,0.50,1.0750,0.333333,0.40,0.0750,1.950,0.9250,0.750,0.20,1.466667,1.150,0.4750,2.60,1.633333,1.30,0.566667,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.30,1.20,1.20,1.0,0.666667,0.050,0.0,0.150,0.50,0.650,0.10,0.850,0.50,0.70,0.0,1.10,1.10,0.40,0.50,0.60,0.60,1.20,1.40,1.20,0.70,2.0,2.50,2.30,1.80,0.60,0.80,0.50,0.348148,0.40,0.80,0.0,0.40,0.50,1.0,0.7750,0.0,0.0,1.0,1.050,1.20,0.80,0.50,1.10,1.10,1.950,0.90,3.40,2.20,0.90,0.60,0.60,0.30,1.10,1.366667,1.10,0.60,1.10,0.70,0.60,1.0,0.70,0.60,0.50,1.0,1.40,0.80,2.266667,1.0,0.966667,0.30,1.40,2.30,2.50,2.10,0.40,0.30,0.50,0.70,1.40,0.50,0.90,1.10,1.10,1.70,0.80,1.0,0.0,2.366667,1.933333,1.066667,2.10,1.40,1.20,0.80,0.60,0.0,0.30,1.40,2.0,2.250,1.30,2.0,1.80,0.50,0.40,1.633333,0.20,0.0,0.0,0.0,0.750,0.30,0.80,0.70,0.60,0.40,2.20,2.50,0.60,1.30,1.50,1.30,1.30,0.455556,0.40,1.40,1.0750,1.70,1.60,0.60,1.20,1.0,0.30,1.60,1.30,1.80,1.20,1.90,0.0,1.866667,1.40,0.90,1.10,1.10,2.40,2.0,1.30,1.60,0.80,0.850,1.30,0.70,1.10,2.30,0.50,1.30,0.30,0.10,0.60,1.10,1.30,2.60,1.70,0.50,0.133333,0.50,0.0,2.0,0.20,0.20,1.80,0.70,1.10,0.30,0.30,0.80,0.30,0.10,1.10,0.90,0.60,2.40,0.60,1.30,2.20,1.10,1.50,1.150,1.40,1.30,0.90,2.0,1.90,1.80,1.50,1.90,1.350,0.50,0.30,1.90,1.550,0.483333,1.20,1.80,2.80,1.30,3.70,2.80,2.40,1.50,0.10,2.30,1.2750,1.80,1.366667,1.30,2.20,0.0,1.80,1.40,1.0,1.30,1.30,0.80,1.40,0.50,0.550,0.50,0.80,0.50,1.40,0.50,0.80,0.70,0.90,0.20,1.0,0.50,0.0,0.10,0.0,0.10,1.0,1.0,0.0,0.60,2.30,0.20,0.10,1.20,1.30,2.10,0.30,0.60,0.70,0.70,0.10,0.840,0.30,0.90,0.90,0.50,0.20,0.20,1.80,2.20,1.20,1.10,1.550,0.80,1.0,0.50,1.70,0.40,1.40,1.20,0.850,0.50,1.30,0.90,0.850,0.60,0.80,1.10,0.0,1.0,1.40,2.20,2.20,1.90,0.0,0.0,1.0,0.30,0.40,1.30,0.60,0.30,1.20,0.90,1.30,1.0,0.20,0.50,0.80,1.10,0.30,1.0,0.60,0.20,1.30,0.60,0.80,1.40,0.80,0.80,2.20,2.60,0.30,2.10,2.30,2.0,0.50,1.50,1.60,1.60,1.20,1.50,0.70,1.0,1.0,1.30,2.90,0.90,0.50,1.30,0.30,1.70,0.90,0.20,0.20,0.30,1.70,0.70,0.40,1.70,1.20,1.70,1.20,0.90,1.70,0.90,1.20,2.10,0.80,2.50,0.30,0.666667,1.0,2.0,1.70,1.60,2.50,2.250,4.10,0.50,0.20,0.10,0.250,0.0,1.450,1.70,1.20,1.50,3.50,3.10,1.10,1.60,1.80,1.350,1.0,2.142857,1.30,1.60,1.750,1.10,1.40,1.20,0.10,1.60,1.60,1.30,0.20,2.20,1.050,0.80,0.40,1.10,0.60,1.50,1.10,0.0,1.585714,1.250,0.70,0.90,1.0,0.0,0.60,1.70,2.750,2.060,1.0,1.70,1.60,2.50,3.30,2.10,1.60,1.20,0.60,1.10,0.90,2.0,0.90,0.70,0.666667,0.70,0.0,1.30,0.60,0.50,1.60,2.20,0.20,0.80,1.20,0.10,2.0,1.80,2.60,1.20,1.60,1.066667,2.20,2.60,0.90,1.40,0.20,0.80,0.0,0.40,0.80,0.30,0.80,0.20,0.60,0.40,1.10,0.883333,1.20,0.0,0.50,0.0,1.30,1.20,0.90,1.233333,1.90,0.80,0.80,1.30,1.50,2.20,0.60,0.90,0.90,0.50,0.20,0.90,0.20,0.90,0.50,0.0,0.60,0.850,0.40,0.30,0.60,0.20,1.30,0.50,0.70,1.0,1.10,0.70,0.80,1.0,0.50,1.80,0.30,2.0,0.90,0.90,0.20,0.80,0.70,1.0,0.30,0.80,1.50,1.50,1.20,1.70,1.30,0.60,1.30,0.30,0.60,0.70,0.10,1.20,0.20,0.90,0.70,0.60,1.80,0.20,0.40,2.30,1.50,0.10,0.0,1.20,1.340,2.0,1.480,1.460,1.633333,0.70,1.020,1.1750,0.30,1.680,1.7750,1.250,0.920,1.060,0.380,1.520,0.860,0.10,1.580,0.50,1.0750,0.333333,0.40,0.0750,1.950,0.9250,0.750,0.20,1.466667,1.150,0.4750,2.60,1.633333,1.30,0.566667,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.70,1.70,0.30,1.20,0.666667,0.050,0.0,0.150,0.50,0.650,0.20,0.850,1.20,0.70,0.0,0.70,1.10,0.90,0.80,1.10,0.90,1.20,1.60,1.10,1.80,2.1250,2.20,1.10,1.70,0.50,1.0,0.20,0.348148,0.40,1.10,0.0,0.50,0.360,0.90,0.60,0.0,0.0,1.0,1.050,1.0,0.50,0.450,0.90,1.2250,1.950,1.2750,3.40,2.20,1.30,1.0,2.60,0.60,1.30,1.366667,0.70,1.40,1.60,1.10,0.60,0.80,1.20,0.40,0.10,0.90,1.30,0.70,2.266667,0.90,0.966667,0.90,1.50,0.50,2.0,1.80,0.40,0.40,0.90,0.90,1.40,0.50,1.10,1.10,1.20,1.250,0.20,1.30,0.0,2.366667,1.933333,1.066667,1.50,1.0,1.20,0.90,0.70,0.20,0.20,1.30,2.20,2.250,1.20,2.0,1.80,1.20,0.357143,1.633333,0.10,0.0,0.0,0.0,0.750,0.30,0.60,0.70,0.60,0.40,2.20,2.50,0.80,1.30,0.666667,1.1250,1.30,0.455556,0.40,1.628571,1.20,1.60,1.30,0.40,1.40,1.0,0.30,2.80,1.30,2.0,1.40,1.90,0.0,1.866667,1.40,0.90,0.80,0.90,1.70,2.0,1.90,1.50,0.70,1.0,0.70,0.70,1.0,0.70,0.30,1.20,0.20,0.30,0.60,1.70,1.90,2.90,1.70,0.50,0.30,0.30,0.0,1.90,0.20,0.20,1.80,0.60,0.60,0.416667,0.30,0.40,0.80,0.10,0.90,1.233333,0.90,1.433333,0.80,1.10,2.50,1.60,1.20,1.150,0.70,1.10,1.40,2.80,1.50,1.0,1.50,1.70,1.350,0.60,0.50,1.385714,1.550,0.60,1.20,0.90,3.30,2.50,2.20,2.70,2.40,1.80,0.333333,2.60,1.2750,1.60,1.70,1.20,2.30,0.10,2.60,1.30,1.20,1.40,1.0,1.30,1.70,0.90,0.550,0.60,0.70,0.20,1.10,0.60,0.80,0.20,0.80,0.20,0.30,0.3250,0.033333,0.60,0.0,0.10,0.40,1.0,0.050,0.60,1.50,0.20,0.10,0.90,1.80,1.60,0.10,1.20,1.0,0.90,0.30,0.840,0.10,1.60,0.90,0.40,0.30,0.20,2.10,1.60,1.50,0.90,1.50,1.433333,0.60,0.40,2.20,0.60,0.80,0.950,0.850,0.314286,1.30,0.866667,0.850,0.40,1.050,1.80,0.0,0.80,1.40,2.20,2.40,1.70,0.0,0.0,0.80,0.10,0.20,1.70,0.6250,0.4250,1.10,1.0,1.30,1.10,0.333333,0.40,1.0,1.30,0.10,0.90,0.60,0.20,0.90,0.10,0.20,1.10,0.90,1.033333,2.20,2.50,0.40,2.30,2.10,2.20,0.50,0.90,1.70,1.60,1.20,1.0,0.70,1.0,0.80,0.80,2.60,0.70,0.50,0.90,0.30,1.40,0.70,0.60,0.60,0.30,1.10,0.70,0.50,1.70,1.414286,1.10,1.30,1.20,1.90,1.10,0.80,2.60,0.80,3.20,0.30,0.90,1.50,2.70,2.30,3.10,1.70,2.250,2.50,0.70,0.10,0.0,0.250,0.0,1.450,1.70,1.40,2.0,4.30,2.80,1.10,1.30,1.80,0.80,1.50,2.10,1.50,2.30,1.750,1.40,2.0,1.40,0.50,1.70,1.40,2.0,1.0,2.50,1.0,0.90,0.533333,0.60,0.50,1.0,1.0,0.10,1.70,1.250,0.633333,0.90,0.20,0.0,1.10,1.528571,2.750,1.70,0.70,1.80,1.60,1.80,2.90,2.30,1.60,2.20,0.50,1.10,1.50,1.70,0.90,1.10,0.666667,0.70,0.0,1.30,0.60,0.50,1.50,1.90,0.10,1.0,2.050,1.60,2.70,2.40,2.10,0.90,1.70,1.40,1.90,2.10,0.80,1.60,1.20,2.10,0.0,0.70,0.80,0.70,0.50,0.366667,1.10,0.40,1.50,0.60,1.40,0.0,0.30,0.0,0.90,1.20,1.10,1.30,1.40,1.0,0.80,0.80,1.30,1.70,0.30,1.20,0.850,0.650,0.30,0.550,0.10,1.50,0.30,0.0,1.50,0.50,0.50,0.30,0.40,0.20,1.50,0.50,0.90,0.80,0.90,0.550,1.150,0.950,0.60,1.916667,0.533333,2.90,0.90,0.60,0.10,0.80,0.90,0.80,0.20,0.80,0.10,1.50,1.50,1.80,1.40,0.60,1.20,0.50,0.80,0.50,0.0,1.70,0.40,0.90,1.40,0.460,1.50,0.0,0.30,2.90,2.90,0.20,0.0,1.30,0.60,1.70,1.10,1.70,1.70,0.70,1.0,0.90,0.30,1.80,1.70,1.30,0.80,0.40,0.40,1.40,0.80,0.10,1.70,0.80,1.0750,0.333333,0.40,0.0750,1.950,0.9250,0.750,0.20,1.466667,1.150,0.4750,2.60,1.633333,1.30,0.566667,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.70,1.40,0.60,1.60,0.666667,0.050,0.0,0.150,0.50,0.650,0.0,0.850,0.70,0.70,0.0,0.60,0.80,0.80,0.60,1.30,0.80,1.20,1.70,1.0,1.30,2.1250,2.40,1.860,2.50,0.90,1.0,0.60,0.348148,0.40,0.70,0.0,0.10,0.360,0.50,0.60,0.0,0.0,1.0,1.050,1.20,1.0,0.450,1.10,1.2250,1.950,1.2750,3.40,2.20,1.50,0.80,3.0,0.557143,1.60,1.366667,0.90,1.40,1.50,1.10,0.30,1.70,1.40,0.0,0.60,0.70,0.90,1.40,2.266667,0.90,0.966667,0.30,1.20,1.720,1.20,1.40,0.40,0.80,0.30,0.90,1.40,0.50,1.028571,1.10,1.30,1.250,0.40,1.0,0.0,2.366667,1.933333,1.066667,1.10,1.60,1.20,0.80,0.533333,0.080,0.180,1.60,2.20,2.250,1.20,2.0,1.80,0.90,0.357143,1.633333,0.10,0.0,0.0,0.0,0.750,0.30,0.7750,0.70,0.60,0.40,2.20,2.50,0.70,1.50,1.313333,1.1250,1.30,0.455556,0.60,1.0,1.0750,1.8250,1.240,0.40,1.16250,1.0,0.30,1.580,1.60,2.0,1.50,1.90,0.0,1.866667,1.30,0.90,1.20,1.028571,2.0,2.0,2.30,1.70,0.70,0.850,1.10,0.70,0.90,1.120,0.50,1.30,0.10,0.0,0.60,1.80,1.40,2.40,2.0,0.40,0.133333,0.30,0.0,2.0,0.20,0.20,1.70,0.60,1.20,0.416667,0.40,0.60,1.0,0.10,1.0,1.233333,0.90,1.433333,0.70,1.10,2.90,1.50,2.10,1.150,0.70,1.40,1.40,2.0,1.60,1.20,1.50,1.40,1.350,0.520,0.70,1.20,1.550,0.50,1.20,1.40,1.80,1.90,2.7750,3.40,2.90,1.90,0.333333,1.50,1.30,1.0,1.366667,1.50,2.10,0.10,4.20,1.40,0.50,2.10,0.70,1.20,1.60,0.60,0.550,0.60,0.783333,0.3250,0.90,0.80,0.80,0.20,0.90,0.20,0.30,0.70,0.033333,0.10,0.20,0.0,0.666667,1.0,0.050,0.60,2.20,0.40,0.10,0.90,1.50,1.90,0.222222,0.70,0.70,0.90,0.40,0.70,0.10,1.60,0.90,0.30,0.30,0.20,1.90,2.20,1.80,1.0,1.550,1.433333,0.80,0.20,3.50,0.4750,1.0,0.30,0.850,0.50,1.30,0.866667,0.850,0.70,1.050,1.10,0.0,0.70,1.40,2.20,1.90,1.20,0.10,0.0,0.70,0.10,0.30,1.0,0.6250,0.4250,1.10,1.0,1.30,1.10,0.333333,0.40,1.50,1.0,0.50,0.90,0.60,0.20,1.40,0.20,0.40,1.70,0.50,1.033333,2.20,2.50,0.30,2.40,1.80,1.90,0.50,1.20,1.70,1.50,1.20,1.0,0.70,1.10,1.0,1.20,2.814286,0.80,0.60,1.10,0.30,1.60,1.30,0.80,0.10,0.30,1.70,0.80,0.50,1.20,1.80,1.70,1.40,1.0,1.30,0.50,0.30,3.50,0.40,3.10,0.30,0.50,1.40,2.40,2.10,2.4750,2.0,2.250,5.20,0.80,0.0,0.10,0.250,0.0,1.450,1.70,1.50,2.70,5.70,2.90,1.10,1.40,1.80,1.70,1.20,2.20,1.50,1.30,1.750,1.10,1.80,0.90,0.10,1.70,1.20,1.40,1.10,2.60,1.30,0.60,0.533333,0.6250,1.0,0.90,0.70,0.10,0.80,1.250,0.633333,0.90,0.20,0.0,1.0,1.40,2.750,2.10,0.60,1.90,2.20,2.30,2.40,2.30,1.60,1.60,0.80,0.90,1.10,1.50,0.90,0.914286,0.666667,0.70,0.0,1.30,0.60,0.90,2.40,1.80,0.166667,0.60,2.050,1.20,2.60,2.40,2.50,1.20,1.90,1.10,2.20,1.90,0.20,1.30,0.780,3.40,0.0,1.0,1.066667,0.70,1.40,0.30,0.60,0.40,1.70,0.80,1.433333,0.30,0.40,0.0,1.50,1.20,0.90,0.90,1.60,0.70,0.80,0.90,1.40,1.90,0.50,2.0,0.850,0.650,0.20,0.550,0.10,1.30,0.30,0.0,0.90,0.850,0.90,0.30,0.50,0.20,1.666667,0.50,0.80,1.30,0.966667,0.550,1.150,0.950,0.70,2.0,0.50,1.866667,0.850,0.750,0.10,1.0,1.0,1.0,0.30,1.0,0.666667,1.50,1.80,1.70,1.10,0.60,1.20,0.70,0.30,0.60,0.10,1.30,0.30,0.90,1.050,0.50,1.40,0.10,0.350,2.90,1.10,0.10,0.0,1.30,0.90,2.10,0.90,1.30,2.0,0.70,1.0,2.20,0.30,1.90,1.7750,0.90,0.80,1.60,0.40,1.40,0.40,0.0,1.20,0.20,0.90,0.50,0.40,0.0,1.80,1.10,0.80,0.20,1.40,1.10,0.60,2.60,2.0,1.30,0.50,1.20,2.60,0.70,1.80,1.350,1.766667,0.70,0.50,0.50,1.0,1.10,0.033333,0.40,0.80,0.60,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.80,1.30,0.90,1.20,0.666667,0.050,0.0,0.150,0.50,0.650,0.0,0.850,0.80,0.70,0.0,0.90,0.70,0.90,0.60,1.40,1.0,1.40,1.60,1.0,0.80,2.1250,1.70,1.860,2.30,0.90,1.0,0.416667,0.348148,0.40,1.30,0.10,0.10,0.360,0.50,1.0,0.0,0.0,1.0,1.050,0.90,0.90,0.450,1.20,1.2250,1.950,1.2750,3.40,2.20,1.40,0.90,1.70,0.60,1.80,1.366667,0.90,1.40,1.20,1.50,1.30,1.20,1.083333,0.50,0.70,1.40,1.40,0.90,2.266667,1.10,0.966667,0.70,2.10,1.720,1.80,2.0,0.40,0.90,0.70,1.30,1.40,0.50,1.20,1.10,1.60,1.250,0.30,1.60,0.0,2.366667,1.933333,1.066667,1.40,1.0,1.20,0.60,0.50,0.080,0.180,1.10,2.40,2.250,1.20,2.0,1.80,0.60,0.30,1.633333,0.166667,0.0,0.0,0.0,0.750,0.30,0.7750,0.70,0.60,0.40,2.20,2.50,0.90,1.0,1.313333,1.1250,1.30,0.455556,0.30,1.628571,1.0750,2.10,1.240,0.40,0.60,1.0,0.30,1.580,1.5750,1.60,1.40,1.90,0.0,1.866667,1.10,0.90,0.90,1.028571,1.10,2.0,1.961905,1.70,0.70,0.850,1.10,0.70,0.60,1.120,0.60,1.40,0.0,0.10,0.60,1.780,1.70,2.60,1.840,0.70,0.133333,0.30,0.0,1.80,0.20,0.20,2.0,0.80,0.40,0.60,0.30,0.70,1.20,0.0,1.50,1.233333,1.10,1.433333,0.60,1.50,2.528571,1.40,1.0,1.150,0.70,1.0,1.20,2.10,1.60,1.50,1.50,1.70,1.350,0.30,0.70,1.50,1.550,0.483333,1.20,1.10,2.10,1.80,2.7750,2.90,2.50,1.80,0.333333,2.40,1.10,1.0,1.366667,1.0,2.10,0.30,2.20,1.0,1.50,1.0,0.90,1.40,2.0,0.80,0.550,0.60,0.80,0.3250,1.10,0.620,0.80,0.40,1.10,0.20,1.0,0.3250,0.033333,0.50,0.20,0.60,0.666667,1.0,0.050,0.60,1.40,0.20,0.0,0.70,0.60,1.90,0.111111,1.0,0.40,0.90,0.60,1.0,0.183333,1.60,0.90,0.60,0.20,0.20,2.10,2.0,1.040,0.30,1.550,1.433333,1.50,0.10,2.20,0.4750,1.0,0.70,0.850,0.10,1.30,0.866667,0.850,0.70,1.050,1.142857,0.0,1.20,1.40,2.20,1.80,1.90,0.0,0.0,0.90,0.30,0.30,1.0,0.50,0.50,0.90,0.80,1.30,1.80,0.333333,0.80,1.10,1.50,0.242857,1.40,0.60,0.20,1.60,0.20,0.40,1.10,0.70,1.033333,2.20,2.40,0.30,2.20,2.0,2.20,0.50,0.80,1.40,2.50,1.20,1.1250,0.70,1.10,1.0,0.90,2.40,0.80,0.50,1.0,0.70,0.90,1.10,0.60,0.40,0.30,1.380,0.70,0.483333,1.30,2.0,1.457143,0.90,1.120,1.40,0.90,0.60,2.30,0.80,3.70,0.30,0.60,1.142857,2.0,1.90,2.4750,2.183333,2.250,3.40,0.614286,0.0,0.0,0.250,0.0,1.450,1.70,1.0,2.60,5.20,2.80,1.10,0.90,1.80,1.350,1.20,2.70,1.50,1.60,1.750,1.0,1.70,1.50,0.260,2.20,1.30,1.40,1.0,2.20,1.050,0.60,0.533333,0.6250,0.820,1.060,0.880,0.10,1.80,1.250,0.40,0.90,0.60,0.0,0.90,1.30,2.750,2.30,0.70,1.80,2.10,2.30,2.90,2.50,1.60,1.70,0.60,0.60,1.0,1.60,0.90,1.0,0.666667,0.70,0.0,1.30,0.60,1.0,2.10,2.60,0.166667,0.80,2.90,1.60,2.90,2.20,1.90,1.0,2.10,1.066667,1.50,2.60,0.70,2.30,0.80,2.10,0.20,0.70,1.066667,1.20,0.70,0.366667,1.0,0.40,2.10,1.20,1.433333,0.10,0.30,0.0,1.30,1.0,0.90,1.30,1.20,0.0,0.80,1.20,1.30,2.10,0.40,1.10,0.850,0.650,0.0,0.550,0.10,1.50,0.60,0.0,1.50,0.850,1.30,0.30,0.50,0.20,1.555556,0.50,0.90,1.40,0.966667,0.550,1.150,0.950,0.70,1.70,0.50,1.866667,0.850,0.750,0.30,1.0,0.90,0.80,0.20,1.40,0.666667,1.50,2.10,1.50,1.50,1.0,1.40,0.70,0.20,0.20,0.0,0.40,0.0,0.90,1.050,0.20,1.50,0.10,0.350,2.50,1.40,0.20,0.0,1.50,1.80,1.90,1.70,1.30,1.20,0.70,1.20,1.1750,0.30,1.80,1.90,1.250,1.0,0.60,0.30,1.70,1.30,0.40,0.90,0.50,1.20,0.40,0.40,0.20,2.10,1.30,0.90,0.20,1.40,1.20,0.30,2.60,1.40,1.30,0.566667,1.20,2.60,0.70,1.80,1.40,1.30,0.60,0.40,0.30,0.60,1.10,0.0,0.60,0.80,0.40,1.90,1.450,1.450,0.40,0.950,1.50,1.30,1.250,0.30,0.650,1.850,1.150,0.0,0.20,1.650,2.0,0.450,1.450,0.40,1.20,0.50,0.950,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.60,1.30,0.757143,1.0,0.666667,0.050,0.0,0.150,0.50,0.650,0.0,0.850,1.40,0.70,0.0,1.0,1.10,0.80,1.0,1.10,0.90,1.30,1.50,1.0,1.20,2.1250,1.90,1.860,2.50,0.90,1.10,0.416667,0.348148,0.40,0.50,0.40,0.60,0.0,0.80,0.90,0.0,0.0,1.0,1.050,1.50,0.80,0.450,0.80,1.2250,1.950,1.2750,3.40,2.20,1.80,0.60,3.0,0.60,1.90,1.366667,1.10,0.90,1.40,1.50,0.70,1.50,1.083333,0.50,0.40,0.90,1.30,0.60,2.266667,0.80,0.966667,0.90,1.31250,1.720,1.66250,1.70,0.40,1.20,1.0,0.70,1.40,0.50,1.028571,1.10,1.33750,1.250,0.51250,1.228571,0.0,2.366667,1.933333,1.066667,1.40,1.10,1.20,0.70,0.533333,0.080,0.180,1.40,1.50,2.60,1.10,2.0,1.80,0.80,0.60,1.633333,0.166667,0.0,0.0,0.0,0.750,0.30,0.7750,0.70,0.60,0.40,2.20,2.50,0.50,1.40,1.313333,1.1250,1.30,0.455556,0.10,1.70,1.0750,1.30,1.240,0.40,1.70,1.0,0.30,1.580,1.70,2.10,1.0,1.90,0.0,1.866667,1.20,0.90,1.50,1.20,1.90,2.0,2.60,1.90,0.70,0.850,1.10,0.70,1.30,1.120,0.30,1.20,0.0,0.10,0.60,1.780,1.70,2.60,1.840,0.716667,0.133333,0.30,0.0,2.0,0.144444,0.20,1.60,0.60,1.10,0.30,0.40,1.10,0.60,0.10,1.0,1.233333,1.10,1.433333,0.60,0.90,2.40,1.50,0.60,1.150,0.60,1.0,1.0,1.70,1.90,1.316667,1.50,1.70,1.350,0.520,1.10,1.30,1.550,0.0,1.20,1.0,2.0,2.20,2.7750,2.933333,2.533333,1.80,0.333333,1.70,1.2750,1.350,1.366667,1.20,2.40,0.083333,2.40,1.414286,0.70,2.0,1.0,1.60,1.30,0.80,0.550,0.60,0.40,0.3250,1.10,0.620,0.80,0.483333,1.057143,0.20,0.50,0.3250,0.033333,0.20,0.10,0.20,0.80,1.0,0.10,0.60,1.90,0.80,0.0,0.90,0.80,1.90,0.138889,0.80,0.80,0.90,0.40,0.840,0.183333,1.60,0.90,0.40,0.40,0.20,2.70,2.50,1.040,0.90,1.550,1.433333,1.30,0.20,2.20,0.4750,0.90,0.80,0.850,0.30,1.30,0.866667,0.850,0.70,1.050,1.20,0.0,0.80,1.40,2.20,2.60,2.10,0.10,0.0,0.828571,0.20,0.250,1.171429,0.6250,0.4250,1.10,0.80,1.30,1.90,0.333333,0.70,1.20,1.40,0.10,1.133333,1.30,0.20,1.50,0.40,0.40,1.50,0.70,1.033333,2.20,2.0,0.30,2.70,1.80,0.90,0.50,0.80,1.50,1.70,1.20,1.1250,0.70,1.10,0.871429,0.70,3.40,1.30,1.0,1.70,0.80,1.50,0.70,0.50,0.40,0.30,1.380,0.70,0.50,1.70,1.90,1.40,1.60,1.120,1.20,0.40,0.816667,2.90,0.80,2.90,0.30,0.30,1.0,2.0,2.30,2.4750,1.70,2.250,4.20,0.40,0.0,0.10,0.250,0.0,1.450,1.70,2.10,2.214286,4.60,2.50,1.10,1.30,1.80,1.350,1.40,2.0,1.50,1.10,1.50,1.70,1.30,1.30,0.260,2.10,1.10,1.80,0.80,2.50,1.050,0.60,0.533333,0.6250,0.820,1.060,0.880,0.10,0.90,1.250,0.70,0.90,0.40,0.0,0.80,1.30,2.750,2.060,0.80,1.50,1.50,2.20,3.40,2.20,1.60,1.30,0.70,1.0,1.40,1.10,0.90,0.90,0.80,0.70,0.0,1.30,0.60,1.10,1.50,2.080,0.166667,0.70,2.80,2.80,2.580,2.10,2.40,0.60,1.40,1.066667,2.60,2.180,0.70,1.50,0.780,2.10,0.0,0.70,1.066667,0.840,0.70,0.366667,1.0,0.40,1.50,0.90,1.433333,0.10,0.340,0.10,0.90,1.0,0.90,1.30,1.40,1.30,0.80,1.10,1.30,1.90,0.50,1.0,0.850,0.650,0.20,0.550,0.10,2.10,0.50,0.0,1.10,0.850,1.20,0.30,0.70,0.20,1.444444,0.50,0.80,0.80,0.966667,0.550,1.150,0.950,0.50,1.90,0.90,1.866667,0.80,0.750,0.10,1.30,1.10,0.80,0.40,1.160,0.666667,1.50,1.740,1.70,1.3250,0.70,1.240,0.60,0.70,0.60,0.0,1.150,0.2250,0.90,1.050,0.40,1.10,0.10,0.350,2.10,1.80,0.20,0.0,1.20,1.50,2.10,1.60,1.30,1.633333,0.70,1.10,0.90,0.30,1.80,1.60,1.30,1.0,1.20,0.30,1.60,0.70,0.0,1.60,0.40,1.50,0.10,0.40,0.0,1.950,0.80,0.70,0.20,1.466667,1.10,0.40,2.60,1.50,1.30,0.60,1.20,2.60,0.70,1.80,1.30,2.0,0.70,0.60,0.70,1.10,1.10,0.10,0.30,0.90,0.80,1.90,1.30,1.60,0.60,1.30,1.50,1.30,1.10,0.30,0.50,2.10,0.80,0.0,0.20,1.50,2.0,0.50,1.50,0.40,1.20,0.60,1.10,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20,0.90,1.30,0.757143,1.0,0.666667,0.050,0.0,0.150,0.50,0.650,0.10,0.850,0.90,0.40,0.0,0.90,0.80,0.60,1.0,1.30,1.0,1.10,1.20,0.90,1.50,2.1250,2.0250,1.860,2.40,1.30,1.30,0.416667,0.348148,0.40,1.20,0.0,0.28750,0.360,0.50,0.40,0.0,0.0,1.0,1.050,1.80,0.70,0.450,0.70,1.2250,1.950,1.2750,3.40,2.20,1.30,0.70,2.80,0.80,2.0,1.366667,1.50,0.70,1.50,1.30,0.80,1.30,1.083333,0.0,0.20,1.40,1.70,0.90,2.266667,1.20,0.966667,0.56250,1.20,1.720,2.0,1.90,0.40,0.542857,1.0,0.60,1.40,0.50,1.0,1.10,1.40,1.250,0.40,1.80,0.0,2.366667,1.933333,1.066667,1.60,1.30,1.20,0.80,0.533333,0.080,0.180,1.30,1.93750,2.0,1.10,2.0,1.80,1.0,0.40,1.633333,0.166667,0.0,0.0,0.0,0.750,0.30,0.7750,0.70,0.60,0.40,2.20,2.50,0.80,1.20,1.313333,1.1250,1.30,0.455556,0.0,2.50,1.0750,2.90,1.240,0.40,1.40,1.0,0.30,1.580,1.60,1.40,1.60,1.90,0.0,1.866667,1.30,0.90,1.30,1.20,1.80,2.0,1.30,1.68750,0.70,0.850,1.10,0.70,0.90,0.40,0.60,1.30,0.0,0.0,0.60,1.780,1.70,2.60,1.840,0.716667,0.133333,0.30,0.0,1.70,0.111111,0.20,1.80,0.40,0.871429,0.60,0.30,0.614286,0.714286,0.10,1.057143,1.233333,0.90,1.433333,0.60,1.40,2.30,1.80,1.20,1.150,0.942857,1.70,1.70,2.50,1.50,1.316667,1.50,1.90,1.350,0.520,0.70,1.30,1.550,0.0,1.20,1.30,2.0,1.40,2.7750,2.933333,2.533333,1.80,0.333333,2.50,1.2750,1.350,1.366667,1.20,1.40,0.083333,2.10,1.40,1.0,1.60,1.20,1.20,1.40,0.80,0.550,0.60,0.783333,0.3250,1.10,0.620,0.80,0.483333,0.70,0.20,0.50,0.3250,0.033333,0.20,0.20,0.20,0.60,1.0,0.050,0.60,1.50,0.357143,0.10,0.871429,0.60,1.90,0.138889,1.40,0.80,0.90,0.30,1.10,0.10,1.60,0.90,0.30,0.60,0.20,1.985714,2.40,1.040,0.50,1.550,1.433333,1.40,0.0,2.20,0.4750,1.10,0.950,0.850,0.40,1.30,0.866667,0.850,0.70,2.0,0.90,0.0,1.0,1.40,2.20,2.20,2.60,0.60,0.0,0.90,0.40,0.250,0.90,0.6250,0.4250,1.10,1.0,1.30,0.70,0.333333,0.642857,1.50,1.40,0.10,1.133333,0.0,0.20,1.442857,0.40,0.50,1.70,0.70,1.033333,2.20,2.70,0.30,2.80,2.20,1.80,0.50,1.60,2.0,2.50,1.20,1.1250,0.70,1.10,1.10,1.10,2.20,1.10,0.90,1.80,0.70,1.90,1.70,0.50,0.60,0.30,1.380,0.70,0.483333,1.70,2.10,1.50,1.70,1.120,1.40,0.90,0.816667,4.0,0.80,2.90,0.30,0.60,1.191837,2.0,1.90,2.4750,2.183333,2.250,3.20,0.20,0.0,0.10,0.250,0.0,1.450,1.70,1.50,3.10,8.20,1.90,1.10,1.40,1.80,1.350,1.60,2.20,1.50,1.60,2.0,1.085714,2.20,1.50,0.260,2.0,1.0,1.40,0.742857,1.80,1.050,0.60,0.533333,0.6250,0.820,1.060,0.880,0.10,1.10,1.250,0.80,0.90,0.50,0.0,1.20,1.50,2.750,2.060,0.70,2.0,2.10,2.60,3.50,2.60,1.60,1.20,0.90,0.80,1.40,1.50,0.90,1.10,0.40,0.70,0.0,1.30,0.60,1.20,1.90,2.080,0.166667,1.30,2.050,3.0,2.580,2.10,2.233333,0.90,2.90,1.066667,1.90,2.180,0.70,1.40,0.90,2.10,0.20,0.70,1.066667,0.840,0.816667,0.366667,0.80,0.40,1.40,0.90,1.433333,0.0,0.340,0.0,0.90,1.80,0.90,1.30,1.466667,0.333333,1.20,0.50,1.30,2.60,0.80,1.166667,0.850,0.650,0.20,0.550,0.10,1.60,0.50,0.0,1.10,0.850,1.30,0.30,0.50,0.20,1.555556,0.50,1.0,1.10,0.966667,0.550,1.150,0.950,0.70,1.90,0.90,0.70,0.80,0.750,0.40,0.980,0.80,0.880,0.20,1.80,0.40,1.50,2.10,1.680,1.3250,0.50,1.10,0.560,0.60,0.70,0.040,1.150,0.2250,0.90,1.050,0.60,1.50,0.10,0.350,2.20,1.30,0.20,0.0,1.30,1.90,2.20,2.10,1.70,1.633333,0.70,0.80,0.70,0.30,1.10,1.90,1.50,1.0,1.50,0.50,1.50,1.10,0.0,2.50,0.60,0.70,0.333333,0.40,0.10,1.950,0.50,0.60,0.20,1.60,1.20,0.60,2.60,1.633333,1.30,0.60,1.20,2.60,0.70,1.80,1.350,2.0,0.80,0.50,0.50,1.30,1.10,0.0,0.30,0.70,0.60,1.90,1.60,1.30,0.20,0.60,1.50,1.30,1.40,0.30,0.80,1.60,1.50,0.0,0.20,1.80,2.0,0.40,1.40,0.40,1.20,0.40,0.80,0.20,1.30,1.0,1.90,0.60,0.70,0.70,0.30,0.40,1.90,0.90,0.90,1.50,1.20};
+
+       // initialize neighborhoods as preprocessed (elsewhere...)
+       int lengthNIix[] = {3,4,4,3,2,3,4,3,4,3,4,3,4,3,3,4,4,3,3,4,4,2,4,4,3,4,3,3,4,3,4,3,4,4,4,4,3,3,4,4,4,3,3,3,3,3,4,3,4,4,3,3,4,4,4,3,2,3,4,4,4,4,3,4,4,4,3,3,3,3,2,2,4,4,3,4,4,3,4,4,4,4,4,4,4,4,4,3,3,3,3,4,3,3,3,4,4,3,4,4,4,4,4,4,2,2,4,4,2,3,4,4,3,4,4,3,4,4,4,4,4,4,3,4,4,3,4,3,3,3,4,3,4,4,4,4,3,4,4,3,3,4,4,3,3,4,3,4,4,3,4,4,4,4,4,4,4,3,4,3,4,3,3,4,4,4,4,4,3,4,4,3,3,4,3,4,3,4,4,2,3,4,2,3,3,3,4,3,4,4,3,3,3,4,4,4,4,4,4,4,3,3,2,3,3,3,4,4,4,4,3,4,3,4,3,4,4,3,3,4,3,4,3,3,3,4,4,3,3,4,4,2,3,3,4,4,2,3,3,4,4,4,4,3,3,3,4,4,4,4,4,4,4,4,4,3,4,3,4,4,4,3,3,3,3,3,3,4,3,4,4,3,3,4,4,4,4,4,4,4,3,4,3,4,3,4,4,4,4,3,3,4,4,4,4,3,3,4,4,4,4,3,3,3,4,4,4,3,4,4,3,3,4,3,4,4,4,4,4,4,4,4,3,4,3,3,4,3,3,3,4,4,3,3,3,3,3,2,4,3,4,4,4,3,2,3,2,3,4,4,3,3,3,4,4,4,4,4,3,4,4,2,3,4,3,4,4,3,3,2,4,3,1,2,1,4,3,4,3,3,4,4,4,4,4,4,4,3,3,4,4,4,4,4,3,4,3,4,3,3,4,4,3,4,4,4,3,3,3,3,3,3,3,3,3,3,4,4,2,2,3,3,3,3,4,4,3,3,2,4,3,4,3,4,4,4,3,4,4,4,3,3,2,2,3,3,3,2,4,4,4,3,3,3,4,3,4,4,3,2,4,3,4,2,3,4,4,4,3,4,3,4,4,3,3,3,4,4,4,4,2,4,4,2,3,4,4,2,4,3,4,2,2,3,3,3,3,4,3,4,4,4,3,3,3,3,4,4,4,4,3,3,2,4,3,3,3,4,4,4,4,4,2,3,4,4,2,2,4,4,4,3,4,3,3,3,4,3,4,4,3,4,4,3,3,3,4,4,4,1,4,3,2,3,3,2,4,3,3,3,4,4,4,3,4,4,4,3,3,3,4,3,3,3,2,4,3,2,3,3,4,4,4,4,3,4,4,4,2,3,4,4,3,4,3,4,3,3,4,4,2,2,4,4,2,4,3,3,3,3,2,4,4,4,4,3,3,3,3,3,4,4,3,3,3};
+       int NIix_tmp[625][4] = 
+       {
+               {79,431,4,-1},
+               {3,568,542,331},
+               {431,4,82,83},
+               {4,1,431,-1},
+               {432,318,-1,-1},
+               {161,40,488,-1},
+               {162,611,488,40},
+               {506,489,373,-1},
+               {9,169,433,240},
+               {433,515,128,-1},
+               {171,244,170,468},
+               {433,515,129,-1},
+               {127,174,242,615},
+               {171,169,245,-1},
+               {174,588,406,-1},
+               {148,20,502,503},
+               {147,20,118,23},
+               {116,187,16,-1},
+               {181,187,178,-1},
+               {21,187,115,184},
+               {21,15,18,512},
+               {184,148,-1,-1},
+               {23,175,419,214},
+               {512,198,22,16},
+               {542,273,201,-1},
+               {150,561,210,49},
+               {420,422,140,-1},
+               {561,49,150,-1},
+               {83,386,504,224},
+               {232,514,234,-1},
+               {514,510,231,333},
+               {234,612,89,-1},
+               {33,238,593,245},
+               {32,8,247,592},
+               {262,260,38,39},
+               {468,243,240,10},
+               {615,243,41,-1},
+               {615,36,41,-1},
+               {615,13,247,464},
+               {615,464,168,163},
+               {466,615,5,261},
+               {263,162,40,-1},
+               {43,447,616,-1},
+               {42,616,48,-1},
+               {303,43,265,-1},
+               {616,265,510,-1},
+               {562,580,440,301},
+               {265,304,48,-1},
+               {44,580,265,303},
+               {416,471,279,26},
+               {578,606,55,-1},
+               {492,270,452,-1},
+               {269,270,492,452},
+               {50,105,596,614},
+               {56,50,280,592},
+               {50,126,107,-1},
+               {54,592,-1,-1},
+               {457,104,577,-1},
+               {64,616,44,65},
+               {62,333,73,469},
+               {62,469,333,290},
+               {69,595,65,300},
+               {335,59,290,-1},
+               {616,70,304,314},
+               {61,303,65,69},
+               {73,131,61,86},
+               {267,544,500,-1},
+               {289,267,137,-1},
+               {303,265,136,-1},
+               {61,303,131,-1},
+               {290,137,-1,-1},
+               {536,613,-1,-1},
+               {290,137,73,304},
+               {136,288,290,138},
+               {555,278,561,-1},
+               {154,519,78,229},
+               {555,279,564,78},
+               {78,429,75,-1},
+               {291,75,472,519},
+               {318,331,525,329},
+               {317,159,81,327},
+               {319,431,317,79},
+               {319,81,2,320},
+               {319,103,326,24},
+               {542,423,504,380},
+               {331,545,319,425},
+               {616,580,333,90},
+               {528,90,237,-1},
+               {543,206,90,-1},
+               {333,31,335,-1},
+               {334,88,62,-1},
+               {344,601,581,618},
+               {618,556,598,-1},
+               {598,618,556,-1},
+               {547,598,97,-1},
+               {96,124,97,94},
+               {97,399,623,496},
+               {350,479,348,-1},
+               {529,481,607,482},
+               {611,101,7,374},
+               {488,372,161,533},
+               {372,488,533,161},
+               {578,606,596,389},
+               {83,590,504,221},
+               {614,274,-1,-1},
+               {596,606,-1,-1},
+               {270,596,606,492},
+               {104,596,566,492},
+               {109,105,-1,-1},
+               {108,492,456,-1},
+               {606,269,596,518},
+               {494,479,349,597},
+               {574,415,502,-1},
+               {112,574,415,502},
+               {112,574,415,502},
+               {574,118,417,-1},
+               {120,411,415,180},
+               {118,502,184,115},
+               {117,182,415,18},
+               {560,116,181,558},
+               {181,116,17,417},
+               {502,178,184,415},
+               {502,184,178,-1},
+               {496,348,497,349},
+               {608,584,623,95},
+               {499,457,515,-1},
+               {127,593,515,567},
+               {459,169,499,-1},
+               {129,169,592,-1},
+               {169,173,238,-1},
+               {139,502,86,64},
+               {501,69,60,-1},
+               {58,315,139,501},
+               {115,73,137,306},
+               {500,520,116,311},
+               {410,131,502,288},
+               {131,410,73,-1},
+               {501,314,520,502},
+               {73,410,574,502},
+               {86,118,135,-1},
+               {416,540,415,-1},
+               {413,470,414,310},
+               {416,140,144,471},
+               {311,116,268,-1},
+               {574,142,575,-1},
+               {574,147,414,119},
+               {558,575,574,-1},
+               {574,16,15,414},
+               {575,147,15,574},
+               {26,144,199,-1},
+               {503,575,210,189},
+               {214,365,416,422},
+               {350,226,472,229},
+               {555,154,430,479},
+               {397,399,555,597},
+               {429,399,564,97},
+               {393,399,398,397},
+               {156,494,479,-1},
+               {432,4,79,555},
+               {432,80,4,-1},
+               {472,82,541,83},
+               {5,263,488,-1},
+               {41,611,7,-1},
+               {506,490,372,37},
+               {165,292,521,594},
+               {164,167,465,521},
+               {292,506,594,366},
+               {164,292,594,465},
+               {508,617,38,-1},
+               {129,13,173,242},
+               {433,242,129,10},
+               {433,10,14,-1},
+               {407,567,498,-1},
+               {174,14,169,406},
+               {173,14,8,-1},
+               {22,197,577,521},
+               {297,282,507,-1},
+               {507,192,297,176},
+               {560,17,118,18},
+               {182,537,-1,-1},
+               {116,500,122,-1},
+               {120,116,18,19},
+               {118,179,-1,-1},
+               {185,116,21,-1},
+               {118,441,21,-1},
+               {410,187,117,-1},
+               {450,444,560,576},
+               {17,259,185,-1},
+               {343,338,75,363},
+               {196,575,364,166},
+               {22,198,418,-1},
+               {199,446,214,-1},
+               {503,559,194,-1},
+               {446,512,503,198},
+               {195,191,199,192},
+               {23,512,198,192},
+               {150,189,293,166},
+               {190,473,589,175},
+               {559,23,473,22},
+               {559,191,151,473},
+               {599,385,542,-1},
+               {599,24,387,-1},
+               {203,562,-1,-1},
+               {202,510,31,-1},
+               {42,265,613,-1},
+               {447,562,595,-1},
+               {562,544,202,42},
+               {528,332,514,510},
+               {214,49,342,149},
+               {49,364,599,150},
+               {150,49,365,-1},
+               {561,142,421,199},
+               {360,415,421,-1},
+               {150,189,212,26},
+               {365,151,49,-1},
+               {361,575,176,503},
+               {450,294,256,422},
+               {590,381,621,-1},
+               {590,621,381,-1},
+               {591,518,542,621},
+               {272,427,542,-1},
+               {223,28,83,386},
+               {553,83,51,-1},
+               {221,389,104,-1},
+               {218,382,275,-1},
+               {624,564,227,497},
+               {624,554,352,495},
+               {354,155,154,-1},
+               {354,227,229,-1},
+               {228,495,291,392},
+               {354,597,305,399},
+               {234,543,-1,-1},
+               {29,462,231,-1},
+               {232,461,544,-1},
+               {232,236,31,461},
+               {237,332,30,29},
+               {232,31,-1,-1},
+               {463,462,31,-1},
+               {593,33,247,-1},
+               {8,593,617,615},
+               {242,170,567,499},
+               {8,247,32,238},
+               {240,433,299,615},
+               {35,567,41,-1},
+               {10,245,169,-1},
+               {244,13,468,-1},
+               {240,40,13,243},
+               {13,615,33,464},
+               {405,506,473,565},
+               {371,619,374,565},
+               {594,619,512,258},
+               {367,594,473,256},
+               {187,257,18,179},
+               {257,178,259,179},
+               {259,258,185,18},
+               {254,259,258,-1},
+               {294,576,445,420},
+               {276,187,487,-1},
+               {254,576,185,250},
+               {254,512,276,474},
+               {615,168,39,567},
+               {467,37,299,-1},
+               {34,464,12,-1},
+               {41,161,467,-1},
+               {595,616,510,-1},
+               {47,449,580,-1},
+               {46,302,267,-1},
+               {266,314,66,562},
+               {470,412,278,-1},
+               {51,516,110,107},
+               {52,106,614,57},
+               {518,590,542,-1},
+               {220,382,275,-1},
+               {274,219,24,389},
+               {273,275,590,389},
+               {606,274,50,224},
+               {257,259,250,487},
+               {291,597,599,305},
+               {597,561,230,574},
+               {519,142,208,76},
+               {50,283,614,-1},
+               {570,577,284,297},
+               {617,176,215,-1},
+               {570,577,287,285},
+               {617,578,570,-1},
+               {570,287,578,521},
+               {617,578,570,297},
+               {578,283,539,238},
+               {73,301,309,265},
+               {66,267,613,-1},
+               {72,63,60,-1},
+               {277,597,472,561},
+               {446,507,473,197},
+               {294,256,189,251},
+               {293,473,190,251},
+               {296,282,299,-1},
+               {589,295,298,-1},
+               {508,176,521,281},
+               {296,617,466,37},
+               {298,617,242,261},
+               {303,265,61,477},
+               {307,288,46,-1},
+               {307,266,313,-1},
+               {310,69,265,-1},
+               {301,308,72,47},
+               {478,471,392,310},
+               {305,312,48,137},
+               {301,574,313,-1},
+               {304,520,113,305},
+               {310,73,304,58},
+               {392,312,306,-1},
+               {302,137,412,-1},
+               {112,574,310,471},
+               {302,520,112,-1},
+               {311,501,70,112},
+               {312,142,306,137},
+               {79,598,527,83},
+               {598,479,432,81},
+               {83,160,568,579},
+               {83,82,2,85},
+               {579,331,569,424},
+               {330,579,569,538},
+               {321,325,423,-1},
+               {569,83,327,505},
+               {423,325,321,-1},
+               {504,83,324,-1},
+               {423,325,424,330},
+               {526,579,525,-1},
+               {526,568,320,-1},
+               {79,526,325,-1},
+               {321,329,423,504},
+               {79,82,320,84},
+               {207,235,544,-1},
+               {613,86,89,-1},
+               {528,90,562,-1},
+               {528,613,62,-1},
+               {491,569,538,-1},
+               {322,534,-1,-1},
+               {343,291,188,201},
+               {341,4,74,-1},
+               {338,291,343,357},
+               {339,343,74,357},
+               {291,343,208,200},
+               {201,359,341,-1},
+               {91,618,-1,-1},
+               {344,530,485,-1},
+               {485,602,-1,-1},
+               {605,94,600,-1},
+               {496,349,479,497},
+               {623,123,460,393},
+               {97,495,623,-1},
+               {97,156,432,-1},
+               {305,354,227,-1},
+               {354,352,227,390},
+               {228,230,392,352},
+               {531,600,584,480},
+               {481,529,530,602},
+               {359,561,192,599},
+               {577,292,49,-1},
+               {357,363,599,192},
+               {570,358,292,49},
+               {486,281,-1,-1},
+               {357,176,281,-1},
+               {486,283,357,362},
+               {366,210,189,-1},
+               {214,199,175,486},
+               {570,617,364,559},
+               {251,368,576,-1},
+               {594,576,369,-1},
+               {166,122,-1,-1},
+               {619,371,248,374},
+               {619,370,374,-1},
+               {100,-1,-1,-1},
+               {374,163,-1,-1},
+               {488,-1,-1,-1},
+               {550,549,545,569},
+               {550,549,538,-1},
+               {538,376,621,220},
+               {516,387,283,-1},
+               {127,589,383,-1},
+               {381,534,79,103},
+               {428,426,377,590},
+               {220,273,542,381},
+               {599,382,387,201},
+               {387,383,343,382},
+               {387,389,382,200},
+               {387,221,599,52},
+               {389,201,386,-1},
+               {606,577,591,-1},
+               {577,273,384,51},
+               {477,353,470,310},
+               {477,315,352,470},
+               {310,354,477,471},
+               {156,494,349,479},
+               {554,438,531,-1},
+               {624,564,608,227},
+               {624,554,495,-1},
+               {154,156,494,97},
+               {97,429,623,-1},
+               {154,155,564,-1},
+               {573,584,608,402},
+               {572,556,608,584},
+               {573,496,608,-1},
+               {508,405,474,594},
+               {405,594,403,473},
+               {248,474,39,297},
+               {585,14,127,-1},
+               {127,587,498,-1},
+               {588,498,129,-1},
+               {588,610,128,-1},
+               {136,500,441,-1},
+               {500,112,122,-1},
+               {574,558,116,-1},
+               {574,471,120,-1},
+               {142,49,147,-1},
+               {113,557,118,-1},
+               {49,140,120,26},
+               {120,113,558,15},
+               {149,145,-1,-1},
+               {146,22,-1,-1},
+               {120,422,26,-1},
+               {199,120,211,-1},
+               {575,214,189,-1},
+               {324,79,103,-1},
+               {84,538,427,331},
+               {621,426,83,545},
+               {427,545,85,-1},
+               {377,426,220,-1},
+               {381,591,-1,-1},
+               {494,154,155,479},
+               {555,399,97,-1},
+               {4,81,317,3},
+               {568,158,555,-1},
+               {171,11,169,174},
+               {522,624,44,352},
+               {436,624,439,475},
+               {624,440,434,-1},
+               {434,43,229,624},
+               {622,478,352,434},
+               {509,435,522,43},
+               {436,46,475,-1},
+               {410,184,502,-1},
+               {410,115,-1,-1},
+               {118,410,-1,-1},
+               {410,116,186,-1},
+               {186,120,444,-1},
+               {199,150,512,-1},
+               {265,88,-1,-1},
+               {511,42,30,332},
+               {42,265,333,612},
+               {216,256,422,445},
+               {591,275,553,-1},
+               {516,614,283,-1},
+               {456,106,458,-1},
+               {578,52,457,617},
+               {513,107,592,-1},
+               {614,459,453,240},
+               {593,127,239,566},
+               {127,240,614,-1},
+               {456,8,-1,-1},
+               {349,230,154,624},
+               {232,237,236,-1},
+               {237,232,543,202},
+               {231,87,-1,-1},
+               {262,240,41,-1},
+               {615,567,163,164},
+               {468,245,261,5},
+               {615,261,263,35},
+               {40,35,161,-1},
+               {47,613,60,73},
+               {392,279,574,-1},
+               {305,49,519,574},
+               {291,154,188,568},
+               {248,197,594,-1},
+               {199,404,594,-1},
+               {522,523,510,-1},
+               {79,432,327,94},
+               {392,44,352,303},
+               {391,42,477,58},
+               {97,429,317,598},
+               {600,482,-1,-1},
+               {480,604,601,581},
+               {581,618,583,582},
+               {530,601,-1,-1},
+               {618,605,529,-1},
+               {605,547,572,581},
+               {363,283,365,214},
+               {259,371,-1,-1},
+               {506,100,533,161},
+               {7,374,615,-1},
+               {374,533,248,99},
+               {336,425,-1,-1},
+               {517,578,-1,-1},
+               {582,509,564,-1},
+               {157,351,74,-1},
+               {229,97,430,-1},
+               {402,584,583,-1},
+               {622,348,482,607},
+               {610,408,457,-1},
+               {593,592,127,240},
+               {134,410,72,116},
+               {137,410,116,131},
+               {121,112,137,-1},
+               {150,192,560,-1},
+               {84,381,326,-1},
+               {551,426,222,-1},
+               {7,488,41,248},
+               {508,292,177,403},
+               {559,539,403,168},
+               {434,493,460,266},
+               {207,203,616,-1},
+               {616,528,203,-1},
+               {23,259,-1,-1},
+               {455,592,107,492},
+               {30,544,233,-1},
+               {240,13,615,-1},
+               {452,492,110,-1},
+               {106,614,516,57},
+               {591,271,606,388},
+               {471,279,75,354},
+               {134,314,136,58},
+               {589,567,297,617},
+               {524,42,-1,-1},
+               {475,42,440,-1},
+               {523,522,42,440},
+               {432,317,327,423},
+               {321,79,-1,-1},
+               {323,316,-1,-1},
+               {207,335,87,511},
+               {581,600,484,604},
+               {607,582,483,345},
+               {623,572,91,-1},
+               {582,607,483,602},
+               {490,374,611,-1},
+               {376,381,321,-1},
+               {621,376,590,-1},
+               {137,179,118,112},
+               {137,122,112,-1},
+               {621,620,85,424},
+               {508,281,589,199},
+               {305,142,112,-1},
+               {432,97,160,568},
+               {84,24,219,1},
+               {207,88,232,-1},
+               {335,206,514,-1},
+               {375,538,551,-1},
+               {375,621,538,320},
+               {346,618,485,92},
+               {281,577,361,175},
+               {424,-1,-1,-1},
+               {376,538,549,336},
+               {545,381,491,-1},
+               {381,491,-1,-1},
+               {621,376,569,-1},
+               {623,396,401,-1},
+               {494,430,-1,-1},
+               {97,92,111,532},
+               {142,268,119,-1},
+               {412,49,119,-1},
+               {508,199,404,-1},
+               {178,415,150,558},
+               {25,359,278,338},
+               {264,42,206,334},
+               {590,103,518,-1},
+               {399,225,354,349},
+               {250,404,249,374},
+               {513,596,52,107},
+               {240,465,521,-1},
+               {432,476,1,-1},
+               {323,376,527,-1},
+               {366,361,577,284},
+               {623,608,124,-1},
+               {608,605,95,-1},
+               {123,97,349,-1},
+               {114,145,-1,-1},
+               {150,26,144,178},
+               {487,254,619,-1},
+               {281,389,-1,-1},
+               {284,102,454,-1},
+               {327,318,325,-1},
+               {86,65,289,48},
+               {529,484,481,582},
+               {532,581,485,607},
+               {485,496,482,347},
+               {346,484,401,-1},
+               {408,406,586,609},
+               {406,609,409,588},
+               {406,409,588,609},
+               {408,586,-1,-1},
+               {521,508,243,-1},
+               {219,451,563,428},
+               {451,220,217,103},
+               {238,513,459,-1},
+               {456,457,454,239},
+               {250,368,403,-1},
+               {265,86,205,58},
+               {516,105,457,-1},
+               {75,291,278,-1},
+               {79,479,94,327},
+               {340,200,387,357},
+               {355,547,-1,-1},
+               {483,547,-1,-1},
+               {618,601,483,344},
+               {91,547,532,355},
+               {481,485,-1,-1},
+               {584,496,401,484},
+               {274,224,52,-1},
+               {623,572,529,-1},
+               {572,531,484,-1},
+               {409,610,498,-1},
+               {498,409,-1,-1},
+               {488,263,162,101},
+               {510,202,332,205},
+               {448,205,264,334},
+               {566,57,454,458},
+               {40,38,41,-1},
+               {86,45,469,-1},
+               {285,578,298,-1},
+               {547,484,91,-1},
+               {371,576,250,-1},
+               {375,549,538,491},
+               {217,425,535,553},
+               {607,349,438,-1},
+               {349,396,608,-1},
+               {396,622,436,-1}
+       };
+       int** NIix = (int**)malloc(nrow*sizeof(int*));
+       for (int i=0; i<nrow; i++)
+       {
+               NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int));
+               for (int j=0; j<lengthNIix[i]; j++)
+                       NIix[i][j] = NIix_tmp[i][j];
+       }
+
+       // main call to LL minimization
+       double pcoef = 1.0, h = 1e-2, eps = 1e-3;
+       int maxit = 1e2;
+       Parameters params = getVarsWithConvexOptim_core(
+               M, lengthNIix, NIix, nrow, ncol, pcoef, h, eps, maxit, S_TRUE, S_TRUE);
+
+       // free neighborhoods parameters arrays
+       for (int i=0; i<nrow; i++)
+               free(NIix[i]);
+       free(NIix);
+
+       // free results memory
+       free(params.theta);
+       for (int i=0; i<nrow; i++)
+               free(params.f[i]);
+       free(params.f);
+}

diff --git a/src/tests/t.dijkstra.c b/src/tests/t.dijkstra.c
new file mode 100644 (file)
index 0000000..62bbd30
--- /dev/null
+++ b/src/tests/t.dijkstra.c
@@ -0,0 +1,67 @@
+#include "tests/helpers.h"
+#include "sources/dijkstra.h"
+
+//custom graph 1 (connected!)
+void test_dijkstra1()
+{
+       int n = 10;
+       double distsIn[100] = 
+       {
+               NAN,4.0,6.0,3.0,NAN,1.0,NAN,NAN,NAN,NAN,
+               4.0,NAN,NAN,5.0,NAN,1.0,3.0,NAN,1.0,NAN,
+               6.0,NAN,NAN,NAN,NAN,7.0,1.0,NAN,NAN,NAN,
+               3.0,5.0,NAN,NAN,4.0,NAN,NAN,NAN,NAN,1.0,
+               NAN,NAN,NAN,4.0,NAN,NAN,NAN,2.0,3.0,NAN,
+               1.0,1.0,7.0,NAN,NAN,NAN,NAN,3.0,NAN,NAN,
+               NAN,3.0,1.0,NAN,NAN,NAN,NAN,NAN,1.0,NAN,
+               NAN,NAN,NAN,NAN,2.0,3.0,NAN,NAN,NAN,1.0,
+               NAN,1.0,NAN,NAN,3.0,NAN,1.0,NAN,NAN,NAN,
+               NAN,NAN,NAN,1.0,NAN,NAN,NAN,1.0,NAN,NAN
+       };
+       double* distsOut;
+
+       distsOut = dijkstra_core(distsIn, 0, 10);
+       //as by-hand computed, distances should be as follow
+       double shouldOutput0[10] = {0.0,2.0,5.0,3.0,6.0,1.0,4.0,4.0,3.0,4.0};
+       ASSERT_TRUE(checkEqualV(shouldOutput0, distsOut, n));
+       free(distsOut);
+
+       distsOut = dijkstra_core(distsIn, 7, 10);
+       //as by-hand computed, distances should be as follow
+       double shouldOutput7[10] = {4.0,4.0,7.0,2.0,2.0,3.0,6.0,0.0,5.0,1.0};
+       ASSERT_TRUE(checkEqualV(shouldOutput7, distsOut, n));
+       free(distsOut);
+}
+
+//custom graph 2 (connected!)
+void test_dijkstra2()
+{
+       int n = 10;
+       // same as graph 1 above, but link between 1 and 5 is now 4.0 instead of 1.0
+       double distsIn[100] = 
+       {
+               NAN,4.0,6.0,3.0,NAN,1.0,NAN,NAN,NAN,NAN,
+               4.0,NAN,NAN,5.0,NAN,4.0,3.0,NAN,1.0,NAN,
+               6.0,NAN,NAN,NAN,NAN,7.0,1.0,NAN,NAN,NAN,
+               3.0,5.0,NAN,NAN,4.0,NAN,NAN,NAN,NAN,1.0,
+               NAN,NAN,NAN,4.0,NAN,NAN,NAN,2.0,3.0,NAN,
+               1.0,4.0,7.0,NAN,NAN,NAN,NAN,3.0,NAN,NAN,
+               NAN,3.0,1.0,NAN,NAN,NAN,NAN,NAN,1.0,NAN,
+               NAN,NAN,NAN,NAN,2.0,3.0,NAN,NAN,NAN,1.0,
+               NAN,1.0,NAN,NAN,3.0,NAN,1.0,NAN,NAN,NAN,
+               NAN,NAN,NAN,1.0,NAN,NAN,NAN,1.0,NAN,NAN
+       };
+       double* distsOut;
+
+       distsOut = dijkstra_core(distsIn, 0, 10);
+       //as by-hand computed, distances should be as follow
+       double shouldOutput0[10] = {0.0,4.0,6.0,3.0,6.0,1.0,6.0,4.0,5.0,4.0};
+       ASSERT_TRUE(checkEqualV(shouldOutput0, distsOut, n));
+       free(distsOut);
+
+       distsOut = dijkstra_core(distsIn, 7, 10);
+       //as by-hand computed, distances should be as follow
+       double shouldOutput7[10] = {4.0,6.0,7.0,2.0,2.0,3.0,6.0,0.0,5.0,1.0};
+       ASSERT_TRUE(checkEqualV(shouldOutput7, distsOut, n));
+       free(distsOut);
+}

diff --git a/src/tests/t.kmeansClustering.c b/src/tests/t.kmeansClustering.c
new file mode 100644 (file)
index 0000000..e44e5e5
--- /dev/null
+++ b/src/tests/t.kmeansClustering.c
@@ -0,0 +1,84 @@
+#include "tests/helpers.h"
+#include "sources/kmeansClustering.h"
+#include <math.h>
+
+//easy data: already clustered, one cluster = 1 vertex in equilateral triangle
+void test_kmeansClustering1()
+{
+       int n=99;
+       double* distances = (double*)malloc(n*n*sizeof(double));
+       for (int i=0; i<n*n; i++)
+               distances[i] = 1.0;
+       int clustCount = 3, clustSize = n/clustCount; //33
+
+       for (int kounter=0; kounter<clustCount; kounter++)
+       {
+               //cluster k: kounter*33...(kounter+1)*33-1
+               for (int i=kounter*clustSize; i<(kounter+1)*clustSize; i++)
+               {
+                       for (int j=kounter*clustSize; j<(kounter+1)*clustSize; j++)
+                               distances[i+n*j] = 0.0; //high-density cluster...
+               }
+       }
+
+       //call to clustering algorithm
+       int* clusters = kmeansWithDistances_core(distances, n, clustCount, 10, 100);
+
+       ASSERT_TRUE(countDistinctValues(clusters, n) == clustCount);
+       ASSERT_TRUE(checkClustersProportions(clusters, n, clustCount, 1e-10));
+       free(distances);
+       free(clusters);
+}
+
+//three isotropic (well separated) gaussian clusters
+void test_kmeansClustering2()
+{
+       // generate 2D data
+       int n=99, d=2;
+       double* M = (double*)malloc(n*d*sizeof(double));
+       int clustCount = 3, clustSize = n/clustCount; //33
+
+       double ctrs[3][2] =
+       {
+               {-3.0,-3.0},
+               {0.0,0.0},
+               {3.0,3.0}
+       };
+
+       srand(time(NULL));
+       for (int kounter=0; kounter<clustCount; kounter++)
+       {
+               //cluster k: kounter*33...(kounter+1)*33-1
+               for (int i=kounter*clustSize; i<(kounter+1)*clustSize; i++)
+               {
+                       double U = (double)rand()/RAND_MAX;
+                       double V = (double)rand()/RAND_MAX;
+                       double fact = sqrt(-2*log(U));
+                       M[i+n*0] = ctrs[kounter][0] + fact * cos(2*M_PI*V);
+                       M[i+n*1] = ctrs[kounter][1] + fact * sin(2*M_PI*V);
+               }
+       }
+
+       // compute distances matrix
+       double* distances = (double*)calloc(n*n,sizeof(double));
+       for (int i=0; i<n; i++)
+       {
+               for (int ii=0; ii<n; ii++)
+               {
+                       double distance = 0.0;
+                       for (int j=0; j<d; j++)
+                               distance += (M[i+n*j] - M[ii+n*j])*(M[i+n*j] - M[ii+n*j]);
+                       distances[i+n*ii] = sqrt(distance);
+               }
+       }
+       free(M); //no need for initial data anymore
+
+       //call to clustering algorithm
+       int* clusters = kmeansWithDistances_core(distances, n, clustCount, 10, 100);
+
+       ASSERT_TRUE(countDistinctValues(clusters, n) == clustCount);
+       //test that each cluster accounts for 1/3 of total data, +/- 10%
+       ASSERT_TRUE(checkClustersProportions(clusters, n, clustCount, 0.1));
+       free(distances);
+       free(clusters);
+}

diff --git a/src/tests/t.neighbors.c b/src/tests/t.neighbors.c
new file mode 100644 (file)
index 0000000..f8ea1c2
--- /dev/null
+++ b/src/tests/t.neighbors.c
@@ -0,0 +1,156 @@
+#include "tests/helpers.h"
+#include "sources/neighbors.h"
+#include <cgds/List.h>
+
+int** list2int(List** L, int n)
+{
+       int** I = (int**)malloc(n*sizeof(int*));
+       for (int i=0; i<n; i++)
+       {
+               int listSize = list_size(L[i]);
+               I[i] = (int*)malloc((1+listSize)*sizeof(int));
+               I[i][0] = listSize;
+               ListIterator* iterJ = list_get_iterator(L[i]);
+               for (int j=1; j<=listSize; j++)
+               {
+                       listI_get(iterJ, I[i][j]);
+                       listI_move_next(iterJ);
+               }
+               listI_destroy(iterJ);
+       }
+       return I;
+}
+
+//10 lines 12 columns, only NA except on antidiagonal
+// ==> distances computed with coordinates only
+void test_neighbors1()
+{
+       int n = 10, m=12;
+       double M[120] = 
+       {
+               NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,1.0,0.0,0.0,
+               NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,1.0,NAN,1.0,1.0,
+               NAN,NAN,NAN,NAN,NAN,NAN,NAN,1.0,NAN,NAN,2.0,2.0,
+               NAN,NAN,NAN,NAN,NAN,NAN,1.0,NAN,NAN,NAN,3.0,3.0,
+               NAN,NAN,NAN,NAN,NAN,1.0,NAN,NAN,NAN,NAN,4.0,4.0,
+               NAN,NAN,NAN,NAN,1.0,NAN,NAN,NAN,NAN,NAN,5.0,5.0,
+               NAN,NAN,NAN,1.0,NAN,NAN,NAN,NAN,NAN,NAN,6.0,6.0,
+               NAN,NAN,1.0,NAN,NAN,NAN,NAN,NAN,NAN,NAN,7.0,7.0,
+               NAN,1.0,NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,8.0,8.0,
+               1.0,NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,NAN,9.0,9.0
+       };
+
+       double alphas[4] = {-1.0, 0.0, 0.5, 1.0};
+       int k = 2; // no need for more
+       for (int j=0; j<4; j++)
+       {
+               double alpha = alphas[j]; // no impact
+               for (int gmode=0; gmode<4; gmode++)
+               {
+                       List** L = getNeighbors_core(M, alpha, k, gmode, S_FALSE, n, m);
+                       int** NIix = list2int(L, n);
+                       for (int jj=0; jj<n; jj++)
+                               list_destroy(L[jj]);
+                       free(L);
+                       for (int jj=1; jj<n-1; jj++)
+                       {
+                               // interior points: 2 neighbors, below and above [except for gmode==1 and jj==2 or n-2]
+                               if (gmode==1 && (jj==2 || jj==n-3))
+                               {
+                                       ASSERT_TRUE(
+                                               NIix[jj][0] == 3 && (NIix[jj][1] == jj-1 || NIix[jj][2] == jj-1 || NIix[jj][3] == jj-1));
+                                       ASSERT_TRUE(
+                                               NIix[jj][0] == 3 && (NIix[jj][1] == jj+1 || NIix[jj][2] == jj+1 || NIix[jj][3] == jj+1));
+                                       if (jj==2)
+                                       {
+                                               ASSERT_TRUE(
+                                                       NIix[jj][0] == 3 && (NIix[jj][1] == jj-2 || NIix[jj][2] == jj-2 || NIix[jj][3] == jj-2));
+                                       }
+                                       else if (jj==n-3)
+                                       {
+                                               ASSERT_TRUE(
+                                                       NIix[jj][0] == 3 && (NIix[jj][1] == jj+2 || NIix[jj][2] == jj+2 || NIix[jj][3] == jj+2));
+                                       }
+                               }
+                               else
+                               {
+                                       // one neighb below, one neighb above
+                                       ASSERT_TRUE(
+                                               NIix[jj][0] == 2 && (NIix[jj][1] == jj-1 || NIix[jj][2] == jj-1));
+                                       ASSERT_TRUE(
+                                               NIix[jj][0] == 2 && (NIix[jj][1] == jj+1 || NIix[jj][2] == jj+1));
+                               }
+                       }
+                       // boundary points in mode 1 (augmented kNN) or 2 (normal kNN) also have 2 neighbors
+                       if (gmode==1 || gmode==2)
+                       {
+                               ASSERT_TRUE(NIix[0][0] == 2 && (NIix[0][1] == 1 || NIix[0][2] == 1));
+                               ASSERT_TRUE(NIix[0][0] == 2 && (NIix[0][1] == 2 || NIix[0][2] == 2));
+                               ASSERT_TRUE(NIix[n-1][0] == 2 && (NIix[n-1][1] == n-2 || NIix[n-1][2] == n-2));
+                               ASSERT_TRUE(NIix[n-1][0] == 2 && (NIix[n-1][1] == n-3 || NIix[n-1][2] == n-3));
+                       }
+                       else
+                       {
+                               // in mutual kNN or in 'by-quadrants' modes, they only have 1 neighbor.
+                               ASSERT_TRUE(NIix[0][0] == 1 && NIix[0][1] == 1);
+                               ASSERT_TRUE(NIix[n-1][0] == 1 && NIix[n-1][1] == n-2);
+                       }
+
+                       for (int i=0; i<n; i++)
+                               free(NIix[i]);
+                       free(NIix);
+               }
+       }
+}
+
+//10 lines, 6 columns, some NA's...
+void test_neighbors2()
+{
+       int n = 6, m=10;
+       double M[60] = 
+       {
+               NAN,1.0,0.0,NAN,0.0,NAN,3.0,NAN,0.0,0.0,
+               1.0,2.0,0.0,NAN,NAN,0.0,2.0,3.0,1.0,1.0,
+               2.0,3.0,0.0,NAN,NAN,1.0,1.0,2.0,2.0,2.0,
+               3.0,NAN,0.0,NAN,1.0,NAN,NAN,1.0,3.0,3.0,
+               2.0,3.0,NAN,3.0,1.0,2.0,1.0,NAN,4.0,4.0,
+               1.0,0.0,NAN,2.0,NAN,1.0,3.0,NAN,5.0,5.0
+       };
+
+       List** L = getNeighbors_core(M, 0.5, 3, 2, S_FALSE, n, m);
+       int** NIix = list2int(L, n);
+       for (int j=0; j<n; j++)
+               list_destroy(L[j]);
+       free(L);
+
+       // check neighbors of first row
+       ASSERT_TRUE(NIix[0][0] == 3);
+       for (int j=1; j<=3; j++)
+       {
+               ASSERT_TRUE(NIix[0][1] == j || 
+                       NIix[0][2] == j || NIix[0][3] == j);
+       }
+       for (int i=0; i<n; i++)
+               free(NIix[i]);
+       free(NIix);
+
+       L = getNeighbors_core(M, -1.0, 3, 2, S_FALSE, n, m);
+       NIix = list2int(L, n);
+       for (int j=0; j<n; j++)
+               list_destroy(L[j]);
+       free(L);
+
+       // check neighbors of fifth row
+       ASSERT_TRUE(NIix[4][0] == 3);
+       for (int j=2; j<=5; j++)
+       {
+               if (j == 4)
+                       continue; //not self-neighbor
+               ASSERT_TRUE(NIix[4][1] == j || 
+                       NIix[4][2] == j || NIix[4][3] == j);
+       }
+
+       for (int i=0; i<n; i++)
+               free(NIix[i]);
+       free(NIix);
+}