R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # 1. vdp.mixt: moodien loytyminen eri dimensiolla, naytemaarilla ja komponenteilla
> #   -> ainakin nopea check
> 
> #######################################################################
> 
> # Generate random data from five Gaussians. 
> # Detect modes with vdp-gm. 
> # Plot data points and detected clusters with variance ellipses
> 
> #######################################################################
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> #source("˜/Rpackages/netresponse/netresponse/R/detect.responses.R")
> #source("˜/Rpackages/netresponse/netresponse/R/internals.R")
> #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
> #dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> #########  Generate DATA #############################################
> 
> # Generate Nc components from normal-inverseGamma prior
> 
> set.seed(12346)
> 
> dd <- 3   # Dimensionality of data
> Nc <- 5   # Number of components
> Ns <- 200 # Number of data points
> sd0 <- 3  # component spread
> rgam.shape = 2 # parameters for Gamma distribution 
> rgam.scale = 2 # parameters for Gamma distribution to define precisions
> 
> 
> # Generate means and variances (covariance diagonals) for the components 
> component.means <- matrix(rnorm(Nc*dd, mean = 0, sd = sd0), nrow = Nc, ncol = dd)
> component.vars <- matrix(1/rgamma(Nc*dd, shape = rgam.shape, scale = rgam.scale), 
+ 	                 nrow = Nc, ncol = dd)
> component.sds <- sqrt(component.vars)
> 
> 
> # Size for each component -> sample randomly for each data point from uniform distr.
> # i.e. cluster assignments
> sample2comp <- sample.int(Nc, Ns, replace = TRUE)
> 
> D <- array(NA, dim = c(Ns, dd))
> for (i in 1:Ns)  {
+     # component identity of this sample
+     ci <- sample2comp[[i]]
+     cm <- component.means[ci,]
+     csd <- component.sds[ci,]
+     D[i,] <- rnorm(dd, mean = cm, sd = csd)
+ }
> 
> 
> ######################################################################
> 
> # Fit mixture model
> out <- mixture.model(D, mixture.method = "bic")
> 
> # FIXME rowmeans(qofz) is constant but not 1
> #qofz <- P.r.s(t(D), list(mu = out$mu, sd = out$sd, w = out$w), log = FALSE)
> 
> ############################################################
> 
> # Compare input data and results
> 
> ord.out <- order(out$mu[,1])
> ord.in <- order(component.means[,1])
> 
> means.out <- out$mu[ord.out,]
> means.in <- component.means[ord.in,]
> 
> # Cluster stds and variances
> sds.out <- out$sd[ord.out,]
> sds.in  <- sqrt(component.vars[ord.in,])
> 
> # -----------------------------------------------------------
> 
> vars.out <- sds.out^2
> vars.in <- sds.in^2
> 
> # Check correspondence between input and output
> if (length(means.in) == length(means.out)) {
+    cm <- cor(as.vector(means.in), as.vector(means.out))
+    csd <- cor(as.vector(sds.in), as.vector(sds.out))
+ }
> 
> # Plot results (assuming 2D)
> 
> ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
> 
> plot(D, pch = 20, main = paste("Cor.means:", round(cm,3), "/ Cor.sds:", round(csd,3)), xlim = ran, ylim = ran) 
> for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
> for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
> 
> ######################################################
> 
> #for (ci in 1:nrow(means.out))  {
> #    points(means.out[ci,1], means.out[ci,2], col = "red", pch = 19)
> #    el <- ellipse(matrix(c(vars.out[ci,1],0,0,vars.out[ci,2]),2), centre = means.out[ci,])
> #    lines(el, col = "red") 						  
> #}
> 
> #for (ci in 1:nrow(means.in))  {
> #    points(means.in[ci,1], means.in[ci,2], col = "blue", pch = 19)
> #    el <- ellipse(matrix(c(vars.in[ci,1],0,0,vars.in[ci,2]),2), centre = means.in[ci,])
> #    lines(el, col = "blue") 						  
> #}
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   3.70    0.29    3.98 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # 1. vdp.mixt: moodien loytyminen eri dimensiolla, naytemaarilla ja komponenteilla
> #   -> ainakin nopea check
> 
> #######################################################################
> 
> # Generate random data from five Gaussians. 
> # Detect modes with vdp-gm. 
> # Plot data points and detected clusters with variance ellipses
> 
> #######################################################################
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> #source("˜/Rpackages/netresponse/netresponse/R/detect.responses.R")
> #source("˜/Rpackages/netresponse/netresponse/R/internals.R")
> #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
> #dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> #########  Generate DATA #############################################
> 
> # Generate Nc components from normal-inverseGamma prior
> 
> set.seed(12346)
> 
> dd <- 3   # Dimensionality of data
> Nc <- 5   # Number of components
> Ns <- 200 # Number of data points
> sd0 <- 3  # component spread
> rgam.shape = 2 # parameters for Gamma distribution 
> rgam.scale = 2 # parameters for Gamma distribution to define precisions
> 
> 
> # Generate means and variances (covariance diagonals) for the components 
> component.means <- matrix(rnorm(Nc*dd, mean = 0, sd = sd0), nrow = Nc, ncol = dd)
> component.vars <- matrix(1/rgamma(Nc*dd, shape = rgam.shape, scale = rgam.scale), 
+ 	                 nrow = Nc, ncol = dd)
> component.sds <- sqrt(component.vars)
> 
> 
> # Size for each component -> sample randomly for each data point from uniform distr.
> # i.e. cluster assignments
> sample2comp <- sample.int(Nc, Ns, replace = TRUE)
> 
> D <- array(NA, dim = c(Ns, dd))
> for (i in 1:Ns)  {
+     # component identity of this sample
+     ci <- sample2comp[[i]]
+     cm <- component.means[ci,]
+     csd <- component.sds[ci,]
+     D[i,] <- rnorm(dd, mean = cm, sd = csd)
+ }
> 
> 
> ######################################################################
> 
> # Fit mixture model
> out <- mixture.model(D, mixture.method = "bic")
> 
> # FIXME rowmeans(qofz) is constant but not 1
> #qofz <- P.r.s(t(D), list(mu = out$mu, sd = out$sd, w = out$w), log = FALSE)
> 
> ############################################################
> 
> # Compare input data and results
> 
> ord.out <- order(out$mu[,1])
> ord.in <- order(component.means[,1])
> 
> means.out <- out$mu[ord.out,]
> means.in <- component.means[ord.in,]
> 
> # Cluster stds and variances
> sds.out <- out$sd[ord.out,]
> sds.in  <- sqrt(component.vars[ord.in,])
> 
> # -----------------------------------------------------------
> 
> vars.out <- sds.out^2
> vars.in <- sds.in^2
> 
> # Check correspondence between input and output
> if (length(means.in) == length(means.out)) {
+    cm <- cor(as.vector(means.in), as.vector(means.out))
+    csd <- cor(as.vector(sds.in), as.vector(sds.out))
+ }
> 
> # Plot results (assuming 2D)
> 
> ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
> 
> plot(D, pch = 20, main = paste("Cor.means:", round(cm,3), "/ Cor.sds:", round(csd,3)), xlim = ran, ylim = ran) 
> for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
> for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
> 
> ######################################################
> 
> #for (ci in 1:nrow(means.out))  {
> #    points(means.out[ci,1], means.out[ci,2], col = "red", pch = 19)
> #    el <- ellipse(matrix(c(vars.out[ci,1],0,0,vars.out[ci,2]),2), centre = means.out[ci,])
> #    lines(el, col = "red") 						  
> #}
> 
> #for (ci in 1:nrow(means.in))  {
> #    points(means.in[ci,1], means.in[ci,2], col = "blue", pch = 19)
> #    el <- ellipse(matrix(c(vars.in[ci,1],0,0,vars.in[ci,2]),2), centre = means.in[ci,])
> #    lines(el, col = "blue") 						  
> #}
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   4.26    0.25    4.50 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Test script for the ICMg method
> 
> # Load the package
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> data(osmo) # Load data
> 
> # Set parameters
> C.boost = 1
> alpha = 10
> beta = 0.01
> B.num = 10
> B.size = 10
> S.num = 10  
> S.size = 10
> C = 24
> pm0 = 0
> V0 = 1               
> V = 0.1
> 
> # Run combined ICMg sampler
> res = ICMg.combined.sampler(osmo$ppi, osmo$exp, C, alpha, beta, pm0, V0, V, B.num, B.size, S.num, S.size, C.boost) 
Sampling ICMg2...

nodes:10250links:1711observations:133components:24alpha:10beta:0.01

Sampling200iterationcs

Burnin iterations:100

I: 0

n(z):428447452416410420431437412403401421461414416400426441439454443428400450

m(z):706974676668666281836971677778688277767766695870

I:10

convL:-0.489524390411784n(z):3155372674684912254495131718319259594247199320266559462233207343484412363

convN:-0.0155872079653155m(z):59127361138526908413437441714850284110770314852599576

I:20

convL:-0.412808687896035n(z):4065442494276172283425221795332241561250255323271510416229169293468405397

convN:-0.00501641884923117m(z):60120321148828798613438471685053284510879344755548678

I:30

convL:-0.384552507076769n(z):4175732054166042632955841989281245481225243293246479437231190349397400407

convN:-0.00264453874944989m(z):53120331158830797913537471664953294611085324656588778

I:40

convL:-0.368285883192839n(z):4155881974146872712635401933284247504227261290189501435195180347387441454

convN:-0.00398326624906045m(z):52120361168830788013538471654952274610886324656588779

I:50

convL:-0.356969601806985n(z):4056601824057123022525601916286231579214248263190410406176206387376488396

convN:-0.00229380200049286m(z):52120371148830788213437481674952274610688324856598279

I:60

convL:-0.354523310963738n(z):4266381764367082712626031913252244625183237284217400432161194407358435388

convN:-0.00225399357807469m(z):52121371158830788013437481674952274610687324856598379

I:70

convL:-0.346186234534221n(z):4156331613747293122596451883239229594200235276187386417161207451388461408

convN:-0.00211520325085406m(z):51121371158830788013637481664952274610488314855598283

I:80

convL:-0.346123180216727n(z):3976941543597683142386581817243233651220239271148376381157230488375445394

convN:-0.00568335018167217m(z):51121381148830788013637481674952274610388314755598284

I:90

convL:-0.356454450609498n(z):3877251493607292712447241843244250697233182253132376406172211468362458374

convN:-0.00337144917512422m(z):51121371158830798013637481664952274610388314855598283

I:100

convL:-0.337028978906674n(z):3787011713577002962226661903246230689236189263156373402169219495366454369

convN:-0.00136866554967014m(z):52121371148830798013537481674952274610388314855588383

Sample iterations:100

I:110

convL:-0.339569855315746n(z):3666691483677522902566711881234223707253192283136338387164229503387442372

convN:-0.00121491064310809m(z):51122371148830798013637481664952274610388314855598283

I:120

convL:-0.339432114828041n(z):3936641483927763052446611876235219656245177303147392390151213486322491364

convN:-0.0027931084967018m(z):52121371148830798013537481674952274610488314855598282

I:130

convL:-0.341052662222956n(z):3866831533947853062517091729225228714250182302166362413155202537344425349

convN:-0.00157376428454109m(z):52121371158830798013538481664952274610389314755588283

I:140

convL:-0.349487578349292n(z):4136771603827513052596991755226224728211192276159397410129225524338449361

convN:-0.00205319543420444m(z):52121371148830798113537481664952274610388314855598283

I:150

convL:-0.338986243708884n(z):4006511593687202922226881751225238740226196288185377402151235545354471366

convN:-0.00622396657382997m(z):51120371148830798113838481654952274610389314755588283

I:160

convL:-0.357043038033143n(z):4076131603537402902366961686220224781258194293178376370155220556360475409

convN:-0.00245899847938083m(z):51120371148830798113938481644953274610288314755588383

I:170

convL:-0.343762395142124n(z):3786191523787362932726781762210218757226187307181381397166218556341476361

convN:-0.00249220495332707m(z):52120381158830798113438481674952274610389314655588283

I:180

convL:-0.339660561937968n(z):3776161623747412922566781727213233760250167341189378378162213564350450379

convN:-0.00395288169295578m(z):52119371188830798113537481644951274610388314855598284

I:190

convL:-0.344765746521564n(z):3896321713817183052387391789190232725215167327190400374155200546314464389

convN:-0.00204343352528668m(z):52120371168930798113537481644952274610388314855598283

I:200

convL:-0.350028293876431n(z):3756311564007342972437241764197237738239142296200373375162190564331485397

convN:-0.00523761882910459m(z):52121371148830788113537481664953274610388314855598283

DONE

> 
> # Compute component membership probabilities for the data points
> res$comp.memb <- ICMg.get.comp.memberships(osmo$ppi, res)
> 
> # Compute (hard) clustering for nodes
> res$clustering <- apply(res$comp.memb, 2, which.max)
> 
> proc.time()
   user  system elapsed 
  11.46    0.23   11.68 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Test script for the ICMg method
> 
> # Load the package
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> data(osmo) # Load data
> 
> # Set parameters
> C.boost = 1
> alpha = 10
> beta = 0.01
> B.num = 10
> B.size = 10
> S.num = 10  
> S.size = 10
> C = 24
> pm0 = 0
> V0 = 1               
> V = 0.1
> 
> # Run combined ICMg sampler
> res = ICMg.combined.sampler(osmo$ppi, osmo$exp, C, alpha, beta, pm0, V0, V, B.num, B.size, S.num, S.size, C.boost) 
Sampling ICMg2...

nodes:10250links:1711observations:133components:24alpha:10beta:0.01

Sampling200iterationcs

Burnin iterations:100

I: 0

n(z):429454408436415455394423415425449408409442426409394424417459419488425427

m(z):666565608169806169547570767663656475758485886877

I:10

convL:-0.476087774602886n(z):3215303264355301831322537436637324295230316280237340426374227425387253231

convN:-0.00753771066558864m(z):376957979914665601069977815079542691439747107363355

I:20

convL:-0.407703107646097n(z):2575864053305531894297578435726405309225249269249391367438172424300166225

convN:-0.00773108577689154m(z):3269569610314463591119478795077542795439648108383358

I:30

convL:-0.366806473608057n(z):2365743813495571980313580451764382296176267284230374350514169416275140192

convN:-0.00702848238366276m(z):3268569610414462591119380795077552794439648108383259

I:40

convL:-0.376470155999394n(z):2275614203395911993320553472780380231189234304265376326482175462222152196

convN:-0.237919421313274m(z):3369579410414462581129380795074562794439648110383258

I:50

convL:-0.354749668195152n(z):2385994243495902025312493435821349251194229286246418276500168468211156212

convN:-0.00231504357621112m(z):3268579410514462591119279785276562795439646110383358

I:60

convL:-0.349420751474565n(z):1955864243406781979336464415849370251194253266261405271510170461191190191

convN:-0.00277625317890419m(z):3268579410413962591169479785276562795439646110383258

I:70

convL:-0.352330836978962n(z):1886773993156112055320414489874363228222242286222389217504186439174216220

convN:-0.00388944803734947m(z):3269579410613962591149380785276562794439546110383358

I:80

convL:-0.341113075284565n(z):2086484173096761966339410482947368243202222273269380213483187431150229198

convN:-0.418643151989129m(z):3468559210413862591159480785177562794439547110383361

I:90

convL:-0.335028897775131n(z):1916964003096871897316429473972351266215224267240377214487181465157234202

convN:-0.0027860540107515m(z):3268579510513862591169580795275562794439546110383257

I:100

convL:-0.341182955477666n(z):1987084143187081905332439454940362249211227270230370202553150440156214200

convN:-0.00429827494263483m(z):3269579710413962591179479795275562794439546109383256

Sample iterations:100

I:110

convL:-0.333957760906659n(z):2006864032927461911322453408941379246222219275229373221535166438159200226

convN:-0.00668668621327485m(z):3269559610514062591169379785076572794439547110383258

I:120

convL:-0.325458108984111n(z):2047454273147791824305488437995378222202243286207341194495148420171210215

convN:-0.00520315417975964m(z):3268559610514062591179479785077572793439547109383258

I:130

convL:-0.325171220561461n(z):2187624112988291785343436437944384235194222289201357190556159410166191233

convN:-0.0104288282623106m(z):3267559610414262591179478785076572793439647109383259

I:140

convL:-0.326060359904518n(z):2157704183108491738341507468890376220208206271214353181569158388176201223

convN:-0.00542622234443698m(z):3268559710513962591169480785076572794439547109383258

I:150

convL:-0.338744325467224n(z):2187294403167991732389453504947386238177213260269370163531164365146204237

convN:-0.216101973124482m(z):3468579410514062591159579795075572694439547110383257

I:160

convL:-0.331441067384899n(z):2217654313118661740321471490931362237196204258260371168560145393140205204

convN:-0.00599392224531788m(z):3268569510514162591169480795075572793439547109383258

I:170

convL:-0.345103100719442n(z):2167884382807931698346481578940365235219212256249372158489158406159210204

convN:-0.00621294979114691m(z):3269569410613962591159479795075582795439546110383258

I:180

convL:-0.338850167553284n(z):2137474373077821729367460511935364254214222257246374157507137438149207236

convN:-0.00485006998269531m(z):3267579410614162601149280795075572795439547110383357

I:190

convL:-0.333291439956426n(z):2447894312917961687316421476977366266207222252244400133525131452150230244

convN:-0.00868538039706697m(z):3268579310513862601159480795175582796439545110383258

I:200

convL:-0.334494177465643n(z):22779442727978416743674145201046360243190223271219386173507166408143203226

convN:-0.00912103103496033m(z):3368579310513958601159279795175602796439448109383359

DONE

> 
> # Compute component membership probabilities for the data points
> res$comp.memb <- ICMg.get.comp.memberships(osmo$ppi, res)
> 
> # Compute (hard) clustering for nodes
> res$clustering <- apply(res$comp.memb, 2, which.max)
> 
> proc.time()
   user  system elapsed 
   9.17    0.28    9.43 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Validate mixture models
> 
> # Generate random data from five Gaussians. 
> # Detect modes 
> # Plot data points and detected clusters 
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> #fs <- list.files("˜/Rpackages/netresponse/netresponse/R/", full.names = TRUE); for (f in fs) {source(f)}; dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> #########  Generate DATA #######################
> 
> res <- generate.toydata()
> D <- res$data
> component.means <- res$means
> component.sds   <- res$sds
> sample2comp     <- res$sample2comp
> 
> ######################################################################
> 
> par(mfrow = c(2,1))
> 
> for (mm in c("vdp", "bic")) {
+ 
+   # Fit nonparametric Gaussian mixture model
+   #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
+   out <- mixture.model(D, mixture.method = mm, max.responses = 10, pca.basis = FALSE)
+ 
+   ############################################################
+ 
+   # Compare input data and results
+ 
+   ord.out <- order(out$mu[,1])
+   ord.in <- order(component.means[,1])
+ 
+   means.out <- out$mu[ord.out,]
+   means.in <- component.means[ord.in,]
+ 
+   # Cluster stds and variances
+   sds.out <- out$sd[ord.out,]
+   vars.out <- sds.out^2
+ 
+   sds.in  <- component.sds[ord.in,]
+   vars.in <- sds.in^2
+ 
+   # Check correspondence between input and output
+   if (length(means.in) == length(means.out)) {
+     cm <- cor(as.vector(means.in), as.vector(means.out))
+     csd <- cor(as.vector(sds.in), as.vector(sds.out))
+   }
+ 
+   # Plot results (assuming 2D)
+   ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
+ 
+   real.modes <- sample2comp
+   obs.modes <- apply(out$qofz, 1, which.max)
+ 
+   # plot(D, pch = 20, main = paste(mm, "/ cor.means:", round(cm,6), "/ Cor.sds:", round(csd,6)), xlim = ran, ylim = ran) 
+   plot(D, pch = real.modes, col = obs.modes, main = paste(mm, "/ cor.means:", round(cm,6), "/ Cor.sds:", round(csd,6)), xlim = ran, ylim = ran) 
+   for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
+   for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
+ 
+ }
> 
> 
> proc.time()
   user  system elapsed 
   2.64    0.15    2.78 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Validate mixture models
> 
> # Generate random data from five Gaussians. 
> # Detect modes 
> # Plot data points and detected clusters 
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> #fs <- list.files("˜/Rpackages/netresponse/netresponse/R/", full.names = TRUE); for (f in fs) {source(f)}; dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> #########  Generate DATA #######################
> 
> res <- generate.toydata()
> D <- res$data
> component.means <- res$means
> component.sds   <- res$sds
> sample2comp     <- res$sample2comp
> 
> ######################################################################
> 
> par(mfrow = c(2,1))
> 
> for (mm in c("vdp", "bic")) {
+ 
+   # Fit nonparametric Gaussian mixture model
+   #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
+   out <- mixture.model(D, mixture.method = mm, max.responses = 10, pca.basis = FALSE)
+ 
+   ############################################################
+ 
+   # Compare input data and results
+ 
+   ord.out <- order(out$mu[,1])
+   ord.in <- order(component.means[,1])
+ 
+   means.out <- out$mu[ord.out,]
+   means.in <- component.means[ord.in,]
+ 
+   # Cluster stds and variances
+   sds.out <- out$sd[ord.out,]
+   vars.out <- sds.out^2
+ 
+   sds.in  <- component.sds[ord.in,]
+   vars.in <- sds.in^2
+ 
+   # Check correspondence between input and output
+   if (length(means.in) == length(means.out)) {
+     cm <- cor(as.vector(means.in), as.vector(means.out))
+     csd <- cor(as.vector(sds.in), as.vector(sds.out))
+   }
+ 
+   # Plot results (assuming 2D)
+   ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
+ 
+   real.modes <- sample2comp
+   obs.modes <- apply(out$qofz, 1, which.max)
+ 
+   # plot(D, pch = 20, main = paste(mm, "/ cor.means:", round(cm,6), "/ Cor.sds:", round(csd,6)), xlim = ran, ylim = ran) 
+   plot(D, pch = real.modes, col = obs.modes, main = paste(mm, "/ cor.means:", round(cm,6), "/ Cor.sds:", round(csd,6)), xlim = ran, ylim = ran) 
+   for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
+   for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
+ 
+ }
> 
> 
> proc.time()
   user  system elapsed 
   3.07    0.14    3.18 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> # Three MODES
> 
> # set.seed(34884)
> set.seed(3488400)
> 
> Ns <- 200
> Nd <- 2
> 
> D3 <- rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+       	    matrix(rnorm(Ns*Nd, mean = 3), ncol = Nd),
+       	    cbind(rnorm(Ns, mean = -3), rnorm(Ns, mean = 3))
+ 	    )
> 
> #X11()
> par(mfrow = c(2,2))
> for (mm in c("vdp", "bic")) {
+   for (pp in c(FALSE, TRUE)) {
+ 
+     # Fit nonparametric Gaussian mixture model
+     out <- mixture.model(D3, mixture.method = mm, pca.basis = pp)
+     plot(D3, col = apply(out$qofz, 1, which.max), main = paste(mm, "/ pca:",  pp)) 
+ 
+   }
+ }
> 
> # VDP is less sensitive than BIC in detecting Gaussian modes (more
> # separation between the clusters needed)
> 
> # pca.basis option is less important for sensitive detection but
> # it will help to avoid overfitting to unimodal features that
> # are not parallel to the axes (unimodal distribution often becomes
> # splitted in two or more clusters in these cases)
> 
> 
> proc.time()
   user  system elapsed 
   7.10    0.18    7.28 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> # Three MODES
> 
> # set.seed(34884)
> set.seed(3488400)
> 
> Ns <- 200
> Nd <- 2
> 
> D3 <- rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+       	    matrix(rnorm(Ns*Nd, mean = 3), ncol = Nd),
+       	    cbind(rnorm(Ns, mean = -3), rnorm(Ns, mean = 3))
+ 	    )
> 
> #X11()
> par(mfrow = c(2,2))
> for (mm in c("vdp", "bic")) {
+   for (pp in c(FALSE, TRUE)) {
+ 
+     # Fit nonparametric Gaussian mixture model
+     out <- mixture.model(D3, mixture.method = mm, pca.basis = pp)
+     plot(D3, col = apply(out$qofz, 1, which.max), main = paste(mm, "/ pca:",  pp)) 
+ 
+   }
+ }
> 
> # VDP is less sensitive than BIC in detecting Gaussian modes (more
> # separation between the clusters needed)
> 
> # pca.basis option is less important for sensitive detection but
> # it will help to avoid overfitting to unimodal features that
> # are not parallel to the axes (unimodal distribution often becomes
> # splitted in two or more clusters in these cases)
> 
> 
> proc.time()
   user  system elapsed 
   4.60    0.21    4.81 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ 
+ library(netresponse)
+ 
+ # SINGLE MODE
+ 
+ # Produce test data that has full covariance
+ # It is expected that
+ # pca.basis = FALSE splits Gaussian with full covariance into two modes
+ # pca.basis = TRUE should detect just a single mode
+ 
+ Ns <- 200
+ Nd <- 2
+ k <- 1.5
+ 
+ D2 <- matrix(rnorm(Ns*Nd), ncol = Nd) %*% rbind(c(1,k), c(k,1))
+ 
+ par(mfrow = c(2,2))
+ for (mm in c("vdp", "bic")) {
+   for (pp in c(FALSE, TRUE)) {
+ 
+     # Fit nonparametric Gaussian mixture model
+     out <- mixture.model(D2, mixture.method = mm, pca.basis = pp)
+     plot(D2, col = apply(out$qofz, 1, which.max), main = paste("mm:" , mm, "/ pp:",  pp)) 
+ 
+   }
+ }
+ 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> proc.time()
   user  system elapsed 
   2.90    0.26    3.15 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ 
+ library(netresponse)
+ 
+ # SINGLE MODE
+ 
+ # Produce test data that has full covariance
+ # It is expected that
+ # pca.basis = FALSE splits Gaussian with full covariance into two modes
+ # pca.basis = TRUE should detect just a single mode
+ 
+ Ns <- 200
+ Nd <- 2
+ k <- 1.5
+ 
+ D2 <- matrix(rnorm(Ns*Nd), ncol = Nd) %*% rbind(c(1,k), c(k,1))
+ 
+ par(mfrow = c(2,2))
+ for (mm in c("vdp", "bic")) {
+   for (pp in c(FALSE, TRUE)) {
+ 
+     # Fit nonparametric Gaussian mixture model
+     out <- mixture.model(D2, mixture.method = mm, pca.basis = pp)
+     plot(D2, col = apply(out$qofz, 1, which.max), main = paste("mm:" , mm, "/ pp:",  pp)) 
+ 
+   }
+ }
+ 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> 
> proc.time()
   user  system elapsed 
   4.15    0.18    4.32 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # Play with different options and check their effect on  running times for bic and vdp 
> 
> skip <- TRUE
> 
> if (!skip) {
+ 
+   Ns <- 100
+   Nd <- 2
+ 
+   set.seed(3488400)
+ 
+   D <- cbind(
+ 
+      	rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+        	      matrix(rnorm(Ns*Nd, mean = 2), ncol = Nd),
+       	      cbind(rnorm(Ns, mean = -1), rnorm(Ns, mean = 3))
+  	    ), 
+ 
+      	rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+        	      matrix(rnorm(Ns*Nd, mean = 2), ncol = Nd),
+       	      cbind(rnorm(Ns, mean = -1), rnorm(Ns, mean = 3))
+  	    )
+ 	    )
+ 
+   rownames(D) <- paste("R", 1:nrow(D), sep = "-")
+   colnames(D) <- paste("C", 1:ncol(D), sep = "-")
+ 
+   ts <- c()
+   for (mm in c("bic", "vdp")) {
+ 
+ 
+     # NOTE: no PCA basis needed with mixture.method = "bic"
+     tt <- system.time(detect.responses(D, verbose = TRUE, max.responses = 5, 
+ 	   		       mixture.method = mm, information.criterion = "BIC", 
+ 			       merging.threshold = 0, bic.threshold = 0, pca.basis = TRUE))
+ 
+     print(paste(mm, ":", round(tt[["elapsed"]], 3)))
+     ts[[mm]] <- tt[["elapsed"]]
+   }
+ 
+    print(paste(names(ts)[[1]], "/", names(ts)[[2]], ": ", round(ts[[1]]/ts[[2]], 3)))
+ 
+ }
> 
> # -> VDP is much faster when sample sizes increase 
> # 1000 samples -> 25-fold speedup with VDP
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.25    0.03    0.28 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # Play with different options and check their effect on  running times for bic and vdp 
> 
> skip <- TRUE
> 
> if (!skip) {
+ 
+   Ns <- 100
+   Nd <- 2
+ 
+   set.seed(3488400)
+ 
+   D <- cbind(
+ 
+      	rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+        	      matrix(rnorm(Ns*Nd, mean = 2), ncol = Nd),
+       	      cbind(rnorm(Ns, mean = -1), rnorm(Ns, mean = 3))
+  	    ), 
+ 
+      	rbind(matrix(rnorm(Ns*Nd, mean = 0), ncol = Nd), 
+        	      matrix(rnorm(Ns*Nd, mean = 2), ncol = Nd),
+       	      cbind(rnorm(Ns, mean = -1), rnorm(Ns, mean = 3))
+  	    )
+ 	    )
+ 
+   rownames(D) <- paste("R", 1:nrow(D), sep = "-")
+   colnames(D) <- paste("C", 1:ncol(D), sep = "-")
+ 
+   ts <- c()
+   for (mm in c("bic", "vdp")) {
+ 
+ 
+     # NOTE: no PCA basis needed with mixture.method = "bic"
+     tt <- system.time(detect.responses(D, verbose = TRUE, max.responses = 5, 
+ 	   		       mixture.method = mm, information.criterion = "BIC", 
+ 			       merging.threshold = 0, bic.threshold = 0, pca.basis = TRUE))
+ 
+     print(paste(mm, ":", round(tt[["elapsed"]], 3)))
+     ts[[mm]] <- tt[["elapsed"]]
+   }
+ 
+    print(paste(names(ts)[[1]], "/", names(ts)[[2]], ": ", round(ts[[1]]/ts[[2]], 3)))
+ 
+ }
> 
> # -> VDP is much faster when sample sizes increase 
> # 1000 samples -> 25-fold speedup with VDP
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.26    0.03    0.28 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Generate Nc components from normal-inverseGamma prior
> 
> set.seed(12346)
> 
> Ns <- 300
> Nd <- 2
> 
> # Isotropic cloud
> D1 <- matrix(rnorm(Ns*Nd), ncol = Nd) 
> 
> # Single diagonal mode
> D2 <- matrix(rnorm(Ns*Nd), ncol = Nd) %*% rbind(c(1,2), c(2,1)) 
> 
> # Two isotropic modes
> D3 <- rbind(matrix(rnorm(Ns/2*Nd), ncol = Nd), matrix(rnorm(Ns/2*Nd, mean = 3), ncol = Nd))
> D <- cbind(D1, D2, D3)
> 
> colnames(D) <- paste("Feature-",  1:ncol(D), sep = "")
> rownames(D) <- paste("Sample-", 1:nrow(D), sep = "")
> 
> 
> proc.time()
   user  system elapsed 
   0.18    0.04    0.21 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Generate Nc components from normal-inverseGamma prior
> 
> set.seed(12346)
> 
> Ns <- 300
> Nd <- 2
> 
> # Isotropic cloud
> D1 <- matrix(rnorm(Ns*Nd), ncol = Nd) 
> 
> # Single diagonal mode
> D2 <- matrix(rnorm(Ns*Nd), ncol = Nd) %*% rbind(c(1,2), c(2,1)) 
> 
> # Two isotropic modes
> D3 <- rbind(matrix(rnorm(Ns/2*Nd), ncol = Nd), matrix(rnorm(Ns/2*Nd, mean = 3), ncol = Nd))
> D <- cbind(D1, D2, D3)
> 
> colnames(D) <- paste("Feature-",  1:ncol(D), sep = "")
> rownames(D) <- paste("Sample-", 1:nrow(D), sep = "")
> 
> 
> proc.time()
   user  system elapsed 
   0.20    0.06    0.25 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ 
+ # 2. netresponse test
+ # test later with varying parameters
+ 
+ # Load the package
+ library(netresponse)
+ #load("../data/toydata.rda")
+ fs <- list.files("../R/", full.names = TRUE); for (f in fs) {source(f)};
+ 
+ data(toydata)
+ 
+ D <- toydata$emat
+ netw <- toydata$netw
+ 
+ # The toy data is random data with 10 features (genes). 
+ # The features 
+ rf <- c(4, 5, 6)
+ #form a subnetwork with coherent responses
+ # with means 
+ r1 <- c(0, 3, 0)
+ r2 <- c(-5, 0, 2)
+ r3 <- c(5, -3, -3)
+ mu.real <- rbind(r1, r2, r3)
+ # real weights
+ w.real <- c(70, 70, 60)/200
+ # and unit variances
+ rv <- 1
+ 
+ # Fit the model
+ #res <- detect.responses(D, netw, verbose = TRUE, mc.cores = 2)
+ #res <- detect.responses(D, netw, verbose = TRUE, max.responses = 4)
+ 
+ res <- detect.responses(D, netw, verbose = TRUE, max.responses = 3, mixture.method = "bic", information.criterion = "BIC", merging.threshold = 1, bic.threshold = 10, pca.basis = FALSE)
+ 
+ print("OK")
+ 
+ # Subnets (each is a list of nodes)
+ subnets <- get.subnets(res)
+ 
+ # the correct subnet is retrieved in subnet number 2:
+ #> subnet[[2]]
+ #[1] "feat4" "feat5" "feat6"
+ 
+ # how about responses
+ # Retrieve model for the subnetwork with lowest cost function value
+ # means, standard devations and weights for the components
+ if (!is.null(subnets)) {
+ m <- get.model.parameters(res, subnet.id = "Subnet-2")
+ 
+ # order retrieved and real response means by the first feature 
+ # (to ensure responses are listed in the same order)
+ # and compare deviation from correct solution
+ ord.obs <- order(m$mu[,1])
+ ord.real <- order(mu.real[,1])
+ 
+ print(paste("Correlation between real and observed responses:", cor(as.vector(m$mu[ord.obs,]), as.vector(mu.real[ord.real,]))))
+ 
+ # all real variances are 1, compare to observed ones
+ print(paste("Maximum deviation from real variances: ", max(abs(rv - range(m$sd))/rv)))
+ 
+ # weights deviate somewhat, this is likely due to relatively small sample size
+ #print("Maximum deviation from real weights: ")
+ #print( (w.real[ord.real] - m$w[ord.obs])/w.real[ord.real])
+ 
+ print("estimated and real mean matrices")
+ print(m$mu[ord.obs,])
+ print(mu.real[ord.real,])
+ 
+ }
+ 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 8
2 / 8
3 / 8
4 / 8
5 / 8
6 / 8
7 / 8
8 / 8
Compute cost for each variable
Computing model for node 1 / 10
Computing model for node 2 / 10
Computing model for node 3 / 10
Computing model for node 4 / 10
Computing model for node 5 / 10
Computing model for node 6 / 10
Computing model for node 7 / 10
Computing model for node 8 / 10
Computing model for node 9 / 10
Computing model for node 10 / 10
independent models done
Computing delta values for edge  1 / 29 

Computing delta values for edge  2 / 29 

Computing delta values for edge  3 / 29 

Computing delta values for edge  4 / 29 

Computing delta values for edge  5 / 29 

Computing delta values for edge  6 / 29 

Computing delta values for edge  7 / 29 

Computing delta values for edge  8 / 29 

Computing delta values for edge  9 / 29 

Computing delta values for edge  10 / 29 

Computing delta values for edge  11 / 29 

Computing delta values for edge  12 / 29 

Computing delta values for edge  13 / 29 

Computing delta values for edge  14 / 29 

Computing delta values for edge  15 / 29 

Computing delta values for edge  16 / 29 

Computing delta values for edge  17 / 29 

Computing delta values for edge  18 / 29 

Computing delta values for edge  19 / 29 

Computing delta values for edge  20 / 29 

Computing delta values for edge  21 / 29 

Computing delta values for edge  22 / 29 

Computing delta values for edge  23 / 29 

Computing delta values for edge  24 / 29 

Computing delta values for edge  25 / 29 

Computing delta values for edge  26 / 29 

Computing delta values for edge  27 / 29 

Computing delta values for edge  28 / 29 

Computing delta values for edge  29 / 29 

Combining groups,  10  group(s) left...

Combining groups,  9  group(s) left...

Combining groups,  8  group(s) left...

Combining groups,  7  group(s) left...

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

[1] "OK"
[1] "Correlation between real and observed responses: 0.999117848017521"
[1] "Maximum deviation from real variances:  0.0391530538149302"
[1] "estimated and real mean matrices"
           [,1]       [,2]       [,3]
[1,] -4.9334982 -0.1575946  2.1613225
[2,] -0.1299285  3.0047767 -0.1841669
[3,]  5.0738471 -2.9334877 -3.2217492
   [,1] [,2] [,3]
r2   -5    0    2
r1    0    3    0
r3    5   -3   -3
> 
> proc.time()
   user  system elapsed 
  41.01    0.23   41.25 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ 
+ # 2. netresponse test
+ # test later with varying parameters
+ 
+ # Load the package
+ library(netresponse)
+ #load("../data/toydata.rda")
+ fs <- list.files("../R/", full.names = TRUE); for (f in fs) {source(f)};
+ 
+ data(toydata)
+ 
+ D <- toydata$emat
+ netw <- toydata$netw
+ 
+ # The toy data is random data with 10 features (genes). 
+ # The features 
+ rf <- c(4, 5, 6)
+ #form a subnetwork with coherent responses
+ # with means 
+ r1 <- c(0, 3, 0)
+ r2 <- c(-5, 0, 2)
+ r3 <- c(5, -3, -3)
+ mu.real <- rbind(r1, r2, r3)
+ # real weights
+ w.real <- c(70, 70, 60)/200
+ # and unit variances
+ rv <- 1
+ 
+ # Fit the model
+ #res <- detect.responses(D, netw, verbose = TRUE, mc.cores = 2)
+ #res <- detect.responses(D, netw, verbose = TRUE, max.responses = 4)
+ 
+ res <- detect.responses(D, netw, verbose = TRUE, max.responses = 3, mixture.method = "bic", information.criterion = "BIC", merging.threshold = 1, bic.threshold = 10, pca.basis = FALSE)
+ 
+ print("OK")
+ 
+ # Subnets (each is a list of nodes)
+ subnets <- get.subnets(res)
+ 
+ # the correct subnet is retrieved in subnet number 2:
+ #> subnet[[2]]
+ #[1] "feat4" "feat5" "feat6"
+ 
+ # how about responses
+ # Retrieve model for the subnetwork with lowest cost function value
+ # means, standard devations and weights for the components
+ if (!is.null(subnets)) {
+ m <- get.model.parameters(res, subnet.id = "Subnet-2")
+ 
+ # order retrieved and real response means by the first feature 
+ # (to ensure responses are listed in the same order)
+ # and compare deviation from correct solution
+ ord.obs <- order(m$mu[,1])
+ ord.real <- order(mu.real[,1])
+ 
+ print(paste("Correlation between real and observed responses:", cor(as.vector(m$mu[ord.obs,]), as.vector(mu.real[ord.real,]))))
+ 
+ # all real variances are 1, compare to observed ones
+ print(paste("Maximum deviation from real variances: ", max(abs(rv - range(m$sd))/rv)))
+ 
+ # weights deviate somewhat, this is likely due to relatively small sample size
+ #print("Maximum deviation from real weights: ")
+ #print( (w.real[ord.real] - m$w[ord.obs])/w.real[ord.real])
+ 
+ print("estimated and real mean matrices")
+ print(m$mu[ord.obs,])
+ print(mu.real[ord.real,])
+ 
+ }
+ 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 8
2 / 8
3 / 8
4 / 8
5 / 8
6 / 8
7 / 8
8 / 8
Compute cost for each variable
Computing model for node 1 / 10
Computing model for node 2 / 10
Computing model for node 3 / 10
Computing model for node 4 / 10
Computing model for node 5 / 10
Computing model for node 6 / 10
Computing model for node 7 / 10
Computing model for node 8 / 10
Computing model for node 9 / 10
Computing model for node 10 / 10
independent models done
Computing delta values for edge  1 / 29 

Computing delta values for edge  2 / 29 

Computing delta values for edge  3 / 29 

Computing delta values for edge  4 / 29 

Computing delta values for edge  5 / 29 

Computing delta values for edge  6 / 29 

Computing delta values for edge  7 / 29 

Computing delta values for edge  8 / 29 

Computing delta values for edge  9 / 29 

Computing delta values for edge  10 / 29 

Computing delta values for edge  11 / 29 

Computing delta values for edge  12 / 29 

Computing delta values for edge  13 / 29 

Computing delta values for edge  14 / 29 

Computing delta values for edge  15 / 29 

Computing delta values for edge  16 / 29 

Computing delta values for edge  17 / 29 

Computing delta values for edge  18 / 29 

Computing delta values for edge  19 / 29 

Computing delta values for edge  20 / 29 

Computing delta values for edge  21 / 29 

Computing delta values for edge  22 / 29 

Computing delta values for edge  23 / 29 

Computing delta values for edge  24 / 29 

Computing delta values for edge  25 / 29 

Computing delta values for edge  26 / 29 

Computing delta values for edge  27 / 29 

Computing delta values for edge  28 / 29 

Computing delta values for edge  29 / 29 

Combining groups,  10  group(s) left...

Combining groups,  9  group(s) left...

Combining groups,  8  group(s) left...

Combining groups,  7  group(s) left...

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

[1] "OK"
[1] "Correlation between real and observed responses: 0.999117848017521"
[1] "Maximum deviation from real variances:  0.0391530538149302"
[1] "estimated and real mean matrices"
           [,1]       [,2]       [,3]
[1,] -4.9334982 -0.1575946  2.1613225
[2,] -0.1299285  3.0047767 -0.1841669
[3,]  5.0738471 -2.9334877 -3.2217492
   [,1] [,2] [,3]
r2   -5    0    2
r1    0    3    0
r3    5   -3   -3
> 
> proc.time()
   user  system elapsed 
  54.21    0.28   54.54 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ # Visualization
+ 
+ library(netresponse)
+ 
+ #fs <- list.files("˜/Rpackages/netresponse/netresponse/R/", full.names = T); for (f in fs) {source(f)}
+ 
+ source("toydata2.R")
+ 
+ # --------------------------------------------------------------------
+ 
+ set.seed(4243)
+ mixture.method <- "bic"
+ 
+ # --------------------------------------------------------------------
+ 
+ res <- detect.responses(D, verbose = TRUE, max.responses = 10, 
+ 	   		       mixture.method = mixture.method, information.criterion = "BIC", 
+ 			       merging.threshold = 1, bic.threshold = 10, pca.basis = FALSE)
+ 
+ res.pca <- detect.responses(D, verbose = TRUE, max.responses = 10, mixture.method = mixture.method, information.criterion = "BIC", merging.threshold = 1, bic.threshold = 10, pca.basis = TRUE)
+ 
+ # --------------------------------------------------------------------
+ 
+ k <- 1
+ 
+ # Incorrect VDP: two modes detected
+ # Correct BIC: single mode detected
+ subnet.id <- names(get.subnets(res))[[k]]
+ 
+ # Correct: single mode detected (VDP & BIC)
+ subnet.id.pca <- names(get.subnets(res.pca))[[k]]
+ 
+ # --------------------------------------------------------------------------------------------------
+ 
+ vis1 <- plot_responses(res, subnet.id, plot_mode = "pca", main = paste("NoPCA; NoDM"))
+ vis2 <- plot_responses(res, subnet.id, plot_mode = "pca", datamatrix = D, main = "NoPCA, DM")
+ vis3 <- plot_responses(res.pca, subnet.id.pca, plot_mode = "pca", main = "PCA, NoDM")
+ vis4 <- plot_responses(res.pca, subnet.id.pca, plot_mode = "pca", datamatrix = D, main = "PCA, DM")
+ 
+ # With original data: VDP overlearns; BIC works; with full covariance data 
+ # With PCA basis: modes detected ok with both VDP and BIC.
+ 
+ # ------------------------------------------------------------------------
+ 
+ # TODO
+ # pca.plot(res, subnet.id)
+ # plot_subnet(res, subnet.id) 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 5
2 / 5
3 / 5
4 / 5
5 / 5
Compute cost for each variable
Computing model for node 1 / 6
Computing model for node 2 / 6
Computing model for node 3 / 6
Computing model for node 4 / 6
Computing model for node 5 / 6
Computing model for node 6 / 6
independent models done
Computing delta values for edge  1 / 15 

Computing delta values for edge  2 / 15 

Computing delta values for edge  3 / 15 

Computing delta values for edge  4 / 15 

Computing delta values for edge  5 / 15 

Computing delta values for edge  6 / 15 

Computing delta values for edge  7 / 15 

Computing delta values for edge  8 / 15 

Computing delta values for edge  9 / 15 

Computing delta values for edge  10 / 15 

Computing delta values for edge  11 / 15 

Computing delta values for edge  12 / 15 

Computing delta values for edge  13 / 15 

Computing delta values for edge  14 / 15 

Computing delta values for edge  15 / 15 

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

Combining groups,  3  group(s) left...

convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 5
2 / 5
3 / 5
4 / 5
5 / 5
Compute cost for each variable
Computing model for node 1 / 6
Computing model for node 2 / 6
Computing model for node 3 / 6
Computing model for node 4 / 6
Computing model for node 5 / 6
Computing model for node 6 / 6
independent models done
Computing delta values for edge  1 / 15 

Computing delta values for edge  2 / 15 

Computing delta values for edge  3 / 15 

Computing delta values for edge  4 / 15 

Computing delta values for edge  5 / 15 

Computing delta values for edge  6 / 15 

Computing delta values for edge  7 / 15 

Computing delta values for edge  8 / 15 

Computing delta values for edge  9 / 15 

Computing delta values for edge  10 / 15 

Computing delta values for edge  11 / 15 

Computing delta values for edge  12 / 15 

Computing delta values for edge  13 / 15 

Computing delta values for edge  14 / 15 

Computing delta values for edge  15 / 15 

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

Combining groups,  3  group(s) left...

Warning messages:
1: In check.network(network, datamatrix, verbose = verbose) :
  No network provided in function call: assuming fully connected nodes.
2: In check.network(network, datamatrix, verbose = verbose) :
  No network provided in function call: assuming fully connected nodes.
> 
> proc.time()
   user  system elapsed 
  35.00    0.34   35.32 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> skip <- FALSE
> 
> if (!skip) {
+ # Visualization
+ 
+ library(netresponse)
+ 
+ #fs <- list.files("˜/Rpackages/netresponse/netresponse/R/", full.names = T); for (f in fs) {source(f)}
+ 
+ source("toydata2.R")
+ 
+ # --------------------------------------------------------------------
+ 
+ set.seed(4243)
+ mixture.method <- "bic"
+ 
+ # --------------------------------------------------------------------
+ 
+ res <- detect.responses(D, verbose = TRUE, max.responses = 10, 
+ 	   		       mixture.method = mixture.method, information.criterion = "BIC", 
+ 			       merging.threshold = 1, bic.threshold = 10, pca.basis = FALSE)
+ 
+ res.pca <- detect.responses(D, verbose = TRUE, max.responses = 10, mixture.method = mixture.method, information.criterion = "BIC", merging.threshold = 1, bic.threshold = 10, pca.basis = TRUE)
+ 
+ # --------------------------------------------------------------------
+ 
+ k <- 1
+ 
+ # Incorrect VDP: two modes detected
+ # Correct BIC: single mode detected
+ subnet.id <- names(get.subnets(res))[[k]]
+ 
+ # Correct: single mode detected (VDP & BIC)
+ subnet.id.pca <- names(get.subnets(res.pca))[[k]]
+ 
+ # --------------------------------------------------------------------------------------------------
+ 
+ vis1 <- plot_responses(res, subnet.id, plot_mode = "pca", main = paste("NoPCA; NoDM"))
+ vis2 <- plot_responses(res, subnet.id, plot_mode = "pca", datamatrix = D, main = "NoPCA, DM")
+ vis3 <- plot_responses(res.pca, subnet.id.pca, plot_mode = "pca", main = "PCA, NoDM")
+ vis4 <- plot_responses(res.pca, subnet.id.pca, plot_mode = "pca", datamatrix = D, main = "PCA, DM")
+ 
+ # With original data: VDP overlearns; BIC works; with full covariance data 
+ # With PCA basis: modes detected ok with both VDP and BIC.
+ 
+ # ------------------------------------------------------------------------
+ 
+ # TODO
+ # pca.plot(res, subnet.id)
+ # plot_subnet(res, subnet.id) 
+ }
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 5
2 / 5
3 / 5
4 / 5
5 / 5
Compute cost for each variable
Computing model for node 1 / 6
Computing model for node 2 / 6
Computing model for node 3 / 6
Computing model for node 4 / 6
Computing model for node 5 / 6
Computing model for node 6 / 6
independent models done
Computing delta values for edge  1 / 15 

Computing delta values for edge  2 / 15 

Computing delta values for edge  3 / 15 

Computing delta values for edge  4 / 15 

Computing delta values for edge  5 / 15 

Computing delta values for edge  6 / 15 

Computing delta values for edge  7 / 15 

Computing delta values for edge  8 / 15 

Computing delta values for edge  9 / 15 

Computing delta values for edge  10 / 15 

Computing delta values for edge  11 / 15 

Computing delta values for edge  12 / 15 

Computing delta values for edge  13 / 15 

Computing delta values for edge  14 / 15 

Computing delta values for edge  15 / 15 

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

Combining groups,  3  group(s) left...

convert the network into edge matrix
removing self-links
matching the features between network and datamatrix
Filter the network to only keep the edges with highest mutual information
1 / 5
2 / 5
3 / 5
4 / 5
5 / 5
Compute cost for each variable
Computing model for node 1 / 6
Computing model for node 2 / 6
Computing model for node 3 / 6
Computing model for node 4 / 6
Computing model for node 5 / 6
Computing model for node 6 / 6
independent models done
Computing delta values for edge  1 / 15 

Computing delta values for edge  2 / 15 

Computing delta values for edge  3 / 15 

Computing delta values for edge  4 / 15 

Computing delta values for edge  5 / 15 

Computing delta values for edge  6 / 15 

Computing delta values for edge  7 / 15 

Computing delta values for edge  8 / 15 

Computing delta values for edge  9 / 15 

Computing delta values for edge  10 / 15 

Computing delta values for edge  11 / 15 

Computing delta values for edge  12 / 15 

Computing delta values for edge  13 / 15 

Computing delta values for edge  14 / 15 

Computing delta values for edge  15 / 15 

Combining groups,  6  group(s) left...

Combining groups,  5  group(s) left...

Combining groups,  4  group(s) left...

Combining groups,  3  group(s) left...

Warning messages:
1: In check.network(network, datamatrix, verbose = verbose) :
  No network provided in function call: assuming fully connected nodes.
2: In check.network(network, datamatrix, verbose = verbose) :
  No network provided in function call: assuming fully connected nodes.
> 
> proc.time()
   user  system elapsed 
  29.53    0.20   29.71 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # 1. vdp.mixt: moodien loytyminen eri dimensiolla, naytemaarilla ja komponenteilla
> #   -> ainakin nopea check
> 
> #######################################################################
> 
> # Generate random data from five Gaussians. 
> # Detect modes with vdp-gm. 
> # Plot data points and detected clusters with variance ellipses
> 
> #######################################################################
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> #source("˜/Rpackages/netresponse/netresponse/R/detect.responses.R")
> #source("˜/Rpackages/netresponse/netresponse/R/internals.R")
> #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
> #dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> 
> #########  Generate DATA #############################################
> 
> res <- generate.toydata()
> D <- res$data
> component.means <- res$means
> component.sds   <- res$sds
> sample2comp     <- res$sample2comp
> 
> ######################################################################
> 
> # Fit nonparametric Gaussian mixture model
> out <- vdp.mixt(D)
> # out <- vdp.mixt(D, c.max = 3) # try with limited number of components -> OK
> 
> ############################################################
> 
> # Compare input data and results
> 
> ord.out <- order(out$posterior$centroids[,1])
> ord.in <- order(component.means[,1])
> 
> means.out <- out$posterior$centroids[ord.out,]
> means.in <- component.means[ord.in,]
> 
> # Cluster stds and variances
> sds.out <- out$posterior$sds[ord.out,]
> sds.in  <- component.sds[ord.in,]
> vars.out <- sds.out^2
> vars.in <- sds.in^2
> 
> # Check correspondence between input and output
> if (length(means.in) == length(means.out)) {
+    cm <- cor(as.vector(means.in), as.vector(means.out))
+    csd <- cor(as.vector(sds.in), as.vector(sds.out))
+ }
> 
> # Plot results (assuming 2D)
> 
> ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
> 
> plot(D, pch = 20, main = paste("Cor.means:", round(cm,3), "/ Cor.sds:", round(csd,3)), xlim = ran, ylim = ran) 
> for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
> for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
> 
> 
> 
> proc.time()
   user  system elapsed 
   2.48    0.21    2.71 

R version 3.5.1 Patched (2018-07-24 r75005) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # 1. vdp.mixt: moodien loytyminen eri dimensiolla, naytemaarilla ja komponenteilla
> #   -> ainakin nopea check
> 
> #######################################################################
> 
> # Generate random data from five Gaussians. 
> # Detect modes with vdp-gm. 
> # Plot data points and detected clusters with variance ellipses
> 
> #######################################################################
> 
> library(netresponse)
Loading required package: Rgraphviz
Loading required package: graph
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, basename, cbind, colMeans, colSums, colnames,
    dirname, do.call, duplicated, eval, evalq, get, grep, grepl,
    intersect, is.unsorted, lapply, lengths, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, rank, rbind,
    rowMeans, rowSums, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min

Loading required package: grid
Loading required package: minet
Loading required package: mclust
Package 'mclust' version 5.4.1
Type 'citation("mclust")' for citing this R package in publications.
Loading required package: reshape2

netresponse (C) 2008-2016 Leo Lahti et al.

https://github.com/antagomir/netresponse
> #source("˜/Rpackages/netresponse/netresponse/R/detect.responses.R")
> #source("˜/Rpackages/netresponse/netresponse/R/internals.R")
> #source("˜/Rpackages/netresponse/netresponse/R/vdp.mixt.R")
> #dyn.load("/home/tuli/Rpackages/netresponse/netresponse/src/netresponse.so")
> 
> 
> #########  Generate DATA #############################################
> 
> res <- generate.toydata()
> D <- res$data
> component.means <- res$means
> component.sds   <- res$sds
> sample2comp     <- res$sample2comp
> 
> ######################################################################
> 
> # Fit nonparametric Gaussian mixture model
> out <- vdp.mixt(D)
> # out <- vdp.mixt(D, c.max = 3) # try with limited number of components -> OK
> 
> ############################################################
> 
> # Compare input data and results
> 
> ord.out <- order(out$posterior$centroids[,1])
> ord.in <- order(component.means[,1])
> 
> means.out <- out$posterior$centroids[ord.out,]
> means.in <- component.means[ord.in,]
> 
> # Cluster stds and variances
> sds.out <- out$posterior$sds[ord.out,]
> sds.in  <- component.sds[ord.in,]
> vars.out <- sds.out^2
> vars.in <- sds.in^2
> 
> # Check correspondence between input and output
> if (length(means.in) == length(means.out)) {
+    cm <- cor(as.vector(means.in), as.vector(means.out))
+    csd <- cor(as.vector(sds.in), as.vector(sds.out))
+ }
> 
> # Plot results (assuming 2D)
> 
> ran <- range(c(as.vector(means.in - 2*vars.in), 
+                as.vector(means.in + 2*vars.in), 
+ 	       as.vector(means.out + 2*vars.out), 
+ 	       as.vector(means.out - 2*vars.out)))
> 
> plot(D, pch = 20, main = paste("Cor.means:", round(cm,3), "/ Cor.sds:", round(csd,3)), xlim = ran, ylim = ran) 
> for (ci in 1:nrow(means.out))  { add.ellipse(centroid = means.out[ci,], covmat = diag(vars.out[ci,]), col = "red") }
> for (ci in 1:nrow(means.in))  { add.ellipse(centroid = means.in[ci,], covmat = diag(vars.in[ci,]), col = "blue") }
> 
> 
> 
> proc.time()
   user  system elapsed 
   2.73    0.15    2.85