From: Benjamin Auder <benjamin@auder> Date: Tue, 9 Dec 2014 00:44:55 +0000 (+0100) Subject: first commit X-Git-Url: https://git.auder.net/variants/Chakart/css/assets/doc/html/Rlogo.svg?a=commitdiff_plain;h=15d1825d71db60d9b1e653c95feadc4d3d48efd3;p=synclust.git first commit --- 15d1825d71db60d9b1e653c95feadc4d3d48efd3 diff --git a/.gitignore b/.gitignore new file mode 100644 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 index 0000000..e69de29 diff --git a/data/example.RData b/data/example.RData new file mode 100644 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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); +}