blindSourceSeparationInANTsR

play around with simple blind source separation examples with images. see results via this link.

library( ANTsR )

## Loading required package: ANTsRCore

## 
## Attaching package: 'ANTsRCore'

## The following object is masked from 'package:stats':
## 
##     var

## The following objects are masked from 'package:base':
## 
##     all, any, apply, max, min, prod, range, sum

library( fastICA )
library( Matrix )
# recreation of "boxes" dataset in daubechies pnas 2009 and calhoun plos one 2013
V <- matrix(rep(1, 100*100), nrow=100)
v1 <- c(11, 40, 21, 70)
v2 <- c(31, 80, 41, 80)
v1.length <- length(V[v1[1]:v1[2], v1[3]:v1[4]])
c1 <- V
c1[ , ] <- rlogis(length(c(V)), location=-1, scale=1)
c1[v1[1]:v1[2], v1[3]:v1[4]] <- rlogis(v1.length, location=2, scale=1)
c2 <- V
c2[ , ] <- rlogis(length(c(V)), location=-1, scale=1)
v2.length <- length(V[v2[1]:v2[2], v2[3]:v2[4]])
c2[v2[1]:v2[2], v2[3]:v2[4]] <- rlogis(v2.length, location=2, scale=1)
###############################################################
X1 <- 0.5*c1 + 0.5*c2
X2 <- 0.3*c1 + 0.7*c2
X3 <- 0.4*c1 + 0.6*c2
X4 <- 0.8*c1 + 0.2*c2
X5 <- 0.25*c1 + 0.45*c2
X6 <- 0.6*c1 + 0.7*c2
ica.mat <- matrix( 
  c( c(X1), c(X2), 
     c(X3), c(X4),
     c(X5), c(X6) ), nrow=6, byrow=T )

Set up dimensionality reduction.

mask = makeImage( c(100,100), 1 )
ilist = matrixToImages( ica.mat, mask )

Let us look at the raw images.

ct = 1
for ( i in ilist ) {
  print( paste( "Image", ct ) )
  plot( i ); ct = ct + 1
}

## [1] "Image 1"

## [1] "Image 2"

## [1] "Image 3"

## [1] "Image 4"

## [1] "Image 5"

## [1] "Image 6"

Run ICA on the matrix.

nc = 2
myica2 <- fastICA( t( ica.mat ), n.comp = nc, verbose = F )
image(matrix(myica2$S[, nc], nrow=100, byrow=F))
image(matrix(myica2$S[, 2], nrow=100, byrow=F))

image(matrix(myica2$S[, 1], nrow=100, byrow=F))

Try robust PCA>

library( rsvd )
myrpca <- rrpca( (ica.mat), k=2, p=10, q=3, svdalg='rsvd', trace=TRUE)

## 
##  Iteration: 1      k = 2       Fro. error = 0.725665032528178
##  Iteration: 2      k = 2       Fro. error = 0.418573776261823
##  Iteration: 3      k = 2       Fro. error = 0.266177256514672
##  Iteration: 4      k = 3       Fro. error = 0.148911643820819
##  Iteration: 5      k = 3       Fro. error = 0.191231056733063
##  Iteration: 6      k = 3       Fro. error = 0.127042350081028
##  Iteration: 7      k = 3       Fro. error = 0.0549955787221079
##  Iteration: 8      k = 3       Fro. error = 0.0314020339832178
##  Iteration: 9      k = 3       Fro. error = 0.0216255837959094
##  Iteration: 10      k = 3       Fro. error = 0.0145852852460658
##  Iteration: 11      k = 3       Fro. error = 0.00885029697429485
##  Iteration: 12      k = 3       Fro. error = 0.00538776115894
##  Iteration: 13      k = 3       Fro. error = 0.00323041871405128
##  Iteration: 14      k = 3       Fro. error = 0.00188375568369
##  Iteration: 15      k = 3       Fro. error = 0.00107420651134587
##  Iteration: 16      k = 3       Fro. error = 0.000640876302896193

image(matrix(myrpca$S[2, ], nrow=100, byrow=F))

image(matrix(myrpca$S[4,], nrow=100, byrow=F))

Now let us initialize eigenanatomy with the ICA solution and use eanat to refine the results.

################
initlist<-list()
for ( i in 1:nc ) {
  initlist[[ i ]] = makeImage( mask, myica2$S[, i]*-1 )
}
sicamat = icawhiten( ica.mat, 4 )
eanatStruct <- sparseDecom( 
  inmatrix   = sicamat, 
  inmask     = mask,
  sparseness = 0.75, # set this explicitly, negative
  # values allow signed pseudo-eigenvectors in results
  nvecs      = nc,
  smooth     = 1.5,  # components will be smooth
  # negative values use edge-preserving methods (WIP)
  cthresh    = 50,  # get large components
  its        = 25,   # optimize for a "long" time
  mycoption  = 0, maxBased = T )   # this controls orthogonality constraints
  # leave as exercise to test other options
#  initializationList = initlist,
#  priorWeight = 10 , maxBased = T )
eanat = eanatStruct$eigenanatomyimages 
emat1 = matrix( eanat[1,], nrow = 100 ) 
print( paste("sparseness1",sum(abs(eanat[1,])>0)/sum(mask==1) ))

## [1] "sparseness1 1"

image( emat1 )

emat2 = matrix( eanat[2,], nrow = 100 )
image( emat2 )

print( paste("sparseness2",sum(abs(eanat[2,])>0)/sum(mask==1) ))

## [1] "sparseness2 1"

Use our deflation-based algorithm which has fewer necessary parameters and appears to be more stable in this example.

sparsenessParam = (range( abs( myica2$S ) ) * 0.002 )[ 2 ] 
eanat<-eanatDef( 
    inmat      = sicamat,
    mask       = mask, 
    smoother   = 1.0,
    priors     = t(myica2$S),
    sparEpsilon = sparsenessParam ) 
emat1 = matrix( eanat[1,], nrow = 100 ) 
print( paste("sparseness1",sum(abs(eanat[1,])>0)/sum(mask==1) ))

## [1] "sparseness1 0.2228"

image( emat1 )

emat2 = matrix( eanat[2,], nrow = 100 )
image( emat2 )

print( paste("sparseness2",sum(abs(eanat[2,])>0)/sum(mask==1) ))

## [1] "sparseness2 0.2272"

Use our deflation-based algorithm with multi-scale SVD-based initialization. Actually, we just use a single scale here.

dtype = "cov"
spatmat = t( imageDomainToSpatialMatrix( mask, mask ) )
smoomat = knnSmoothingMatrix( spatmat, k = 3^3, sigma = 125 ) # sigma = 250

## Note: method with signature 'CsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
##  target signature 'dsCMatrix#lsCMatrix#missing#numeric'.
##  "Matrix#lsparseMatrix#missing#replValue" would also be valid

icaCorSparse = sparseDistanceMatrix( ica.mat %*% smoomat, k = 1500, kmetric = dtype, eps = 0.0 )
# mysvdFull = irlba::partial_eigen( icaCorSparse, nc*2 ); mysvd = mysvdFull$vec
mysvdFull = irlba::irlba( icaCorSparse, nc*2 )
mysvd = mysvdFull$v
sparsenessParam = ( range( abs( mysvd ) ) * 0.002 )[ 2 ]
sparsenessParam = 1e-4
# how does the svd result - alone - look ?
plot( mask, makeImage( mask , abs( mysvdFull$v[,1] ) ) %>% iMath("Normalize"))

## NULL

plot( mask, makeImage( mask , abs( mysvdFull$v[,2] ) ) %>% iMath("Normalize"))

## NULL

if ( T )
  eanat<-eanatDef( 
      inmat      = scale( ica.mat ),
      mask       = mask, 
      smoother   = 1.0, cthresh=1000,
      priors     = t(mysvd),
      positivity = TRUE,
      sparEpsilon = sparsenessParam ) else eanat = t( mysvd )
for ( i in 1:nrow( eanat ) )
  {
  emat1 = matrix( eanat[i,], nrow = 100 ) 
  print( paste("sparseness1",sum(abs(eanat[i,])>sparsenessParam)/sum(mask==1) ))
  image( emat1 )
  Sys.sleep( 1 )
  }

## [1] "sparseness1 0.205"

## [1] "sparseness1 0.2168"

## [1] "sparseness1 0.8599"

## [1] "sparseness1 0.7894"

Try some neuroimaging data.

img = antsImageRead( getANTsRData( 'r16') )
ilist = list(
  antsImageRead( getANTsRData( 'r27') ),
  antsImageRead( getANTsRData( 'r30') ),
  antsImageRead( getANTsRData( 'r62') ),
  antsImageRead( getANTsRData( 'r64') ),
  antsImageRead( getANTsRData( 'r85') ) )
rlist = list( )
for ( i in 1:length( ilist ) ) {
  temp = antsRegistration( img, ilist[[i]], typeofTransform = "SyN" )
  jacobian = createJacobianDeterminantImage( img, temp$fwdtransforms[[1]], 1)
  rlist[[ i ]] = jacobian 
  rlist[[ i ]] = temp$warpedmov 
  }
imask = getMask( img )
mat   = imageListToMatrix( rlist, imask )

SMP!

spatmat = t( imageDomainToSpatialMatrix( imask, imask ) )
smoomat = knnSmoothingMatrix( spatmat, k = 5^3, sigma = 10 )
myica2 <- fastICA( t( mat ), n.comp = 4, verbose = F )
bdf = data.frame( t( myica2$A ) )
bdf = data.frame( mat %*% t( rsvd::rrpca( scale( mat, scale=F ), nu=4  )$S ) )
myform = as.formula( " x ~ X1 + X2 + X3 + X4  " )
x = scale( mat )
smp = smoothMatrixPrediction( 
  x=x, basisDf=bdf, modelFormula = myform, gamma = 0.01,
  iterations = 25,
  sparsenessQuantile = 0.6, smoothingMatrix = smoomat, 
  positivity = T, verbose = F )

## reducing gradient step: 0.005

## reducing gradient step: 0.0025

## reducing gradient step: 0.00125

## reducing gradient step: 0.000625

par( mfrow = c(1, ncol(smp$v)) )
for ( k in 1:ncol(smp$v) )
  plot( img, makeImage( imask, abs(smp$v[,k]) ) %>% iMath("Normalize"), window.overlay=c(0.1,1) )

par( mfrow = c(1, 1) )
temp = mat %*% smp$v
lmlist = summary( lm( temp  ~ X1 + X2 , data=bdf  ) )
pander::pander( lmlist[[1]] )

	Estimate	Std. Error	t value	Pr(>\|t\|)
X1	-0.0008879	0.0002602	-3.413	0.07619
X2	-0.0001286	0.0002497	-0.515	0.6578
(Intercept)	56111	4711	11.91	0.006976

Fitting linear model: X1 ~ X1 + X2
Observations	Residual Std. Error	\(R^2\)	Adjusted \(R^2\)
5	1618	0.8581	0.7162

pander::pander( lmlist[[2]] )

	Estimate	Std. Error	t value	Pr(>\|t\|)
X1	-0.0002347	0.0001043	-2.25	0.1534
X2	-0.0003575	0.0001001	-3.57	0.07029
(Intercept)	25651	1889	13.58	0.00538

Fitting linear model: X2 ~ X1 + X2
Observations	Residual Std. Error	\(R^2\)	Adjusted \(R^2\)
5	648.8	0.8749	0.7498

Run a population eigenanatomy with multi-res SVD initialization.

# here is an efficient way to make a sparse matrix and get its SVD
myk = 200
icaCorSparse = sparseDistanceMatrix( icawhiten( mat, 4 ) , myk, kmetric=dtype )
mysvdFull = irlba::partial_eigen( icaCorSparse, 4 )
mysvd = mysvdFull$vec
sparsenessParam = ( range( abs( mysvd ) ) * 0.00001 )[ 2 ]
sparsenessParam = 1e-6
if ( T )
  eanat<-eanatDef( 
      inmat      = mat,
      mask       = imask, 
      smoother   = 1.0, its=5,
      priors     = t(mysvd), priorWeight = 0.0,
      sparEpsilon = sparsenessParam ) else {
        eanat = t( mysvd )
        bica = fastICA( t(mat), 4 )
#        eanat = t( bica$S ) 
      }
eilist = matrixToImages( eanat, imask )
eseg = eigSeg( imask, eilist )
plot( img, eseg )

Now try robust pca again …

library( rsvd )
myrpca <- rrpca( mat, k=2, p=10, q=3, svdalg='rsvd', trace=TRUE)

## 
##  Iteration: 1      k = 2       Fro. error = 0.754635934177966
##  Iteration: 2      k = 2       Fro. error = 0.681482876620957
##  Iteration: 3      k = 2       Fro. error = 0.358807863811469
##  Iteration: 4      k = 2       Fro. error = 0.266082228281336
##  Iteration: 5      k = 2       Fro. error = 0.163775136381329
##  Iteration: 6      k = 2       Fro. error = 0.0811873001736234
##  Iteration: 7      k = 2       Fro. error = 0.0426962339002268
##  Iteration: 8      k = 2       Fro. error = 0.0210254779762081
##  Iteration: 9      k = 2       Fro. error = 0.01114436838704
##  Iteration: 10      k = 2       Fro. error = 0.00556569594363045
##  Iteration: 11      k = 2       Fro. error = 0.00271123692248189
##  Iteration: 12      k = 2       Fro. error = 0.00148671054562867
##  Iteration: 13      k = 2       Fro. error = 0.000894938095814845

reilist = matrixToImages( myrpca$S, imask )
eseg = eigSeg( imask, reilist )
k=5
plot( img, reilist[[k]], window.overlay=range(  reilist[[k]] )  )

## NULL

# plot( makeImage(matrix(myrpca$S[2, ], nrow=100, byrow=F))
# image(matrix(myrpca$S[4,], nrow=100, byrow=F))

Try BOLD functional data.

boldfn = getANTsRData( "rsbold" )
bold   = antsImageRead( boldfn )
sbold  = smoothImage( bold, c( rep(6,3),1),
  sigmaInPhysicalCoordinates = TRUE, FWHM = TRUE )
boldavg = getAverageOfTimeSeries( bold )
bmsk   = getMask( boldavg )
plot( boldavg, axis=3 )

Now set up dim red.

set.seed( 11 )
bmat = timeseries2matrix( sbold, bmsk )
cc   = compcor( bmat, 10, variance_extreme = 0.1 )
rbmat = residuals( lm( bmat ~ cc ) )
nc   = 10
# bica = fastICA( t( rbmat ), nc )
myk = 1000
icaCorSparse = sparseDistanceMatrix( rbmat, myk, 
  kmetric=dtype, eps=0.01 )
mysvdFull = irlba::irlba( icaCorSparse, 10 )
bica = mysvdFull$v

View a component.

k = 10
bicaVec = bica[,k]
bicaVec = bicaVec * sign( quantile( bicaVec )[3] )
qbv     = quantile( bicaVec,  0.99 )
bicaVec[ bicaVec > qbv ] = qbv
bicaImg = makeImage( bmsk, bicaVec ) / qbv
plot( boldavg, bicaImg , window.overlay=c(0.1,1), axis=3 )

Eanat with this data.

sparsenessParam = 5.e-3
eanat<-eanatDef( 
    inmat      = scale( rbmat ),
    mask       = bmsk,
    smoother   = 1.0,
    cthresh    = 10,
    positivity = FALSE,
    priors     = t(bica), verbose=T, priorWeight = 0.0,
    sparEpsilon = sparsenessParam ) 
eseg = eigSeg( mask = bmsk, imgList = matrixToImages( eanat, bmsk ) )

View all of the Eanat components in BOLD space.

for ( k in 1:nc )
  {
  tempvec = eanat[k,]
  eimg = makeImage( bmsk, tempvec) / max( tempvec )
  print( k )
  plot( boldavg, eimg , window.overlay=c(0.1,1), axis=3 )
  Sys.sleep( 2 )
  }

Let us look the eigenvalue spectrum over scale.

bmat = timeseries2matrix( bold, bmsk )
tbmat = t( bmat )
bknn  = get.knn( t(bmat), k=100, algo="CR" )
# how do we compute "fast" correlations ...
# 
temp = svd( bmat )$d
dscales = 1:12
dmat = matrix( nrow=length( dscales )+1, 
               ncol = length( temp ) )
dmat[1,] = temp
for ( s in dscales )
  {
  sigma = as.numeric( s )
  sbold  = smoothImage( bold, c( rep(sigma,3), 1.0 ),
    sigmaInPhysicalCoordinates = TRUE, FWHM = TRUE )
  bmat = timeseries2matrix( sbold, bmsk )
  rbmat = residuals( lm( bmat ~ cc ) )
  dmat[ s+1, ] = svd( rbmat )$d
  }
row.names( dmat ) = c( 0, dscales )
pheatmap::pheatmap( dmat[2:(length(dscales)+1),2:10],
                    cluster_rows=F, cluster_cols=F )

The eigenvalues fall off rapidly after scale 3.

Let us investigate the N-sphere. Plot its singular values.

sphereDim = 9
embeddDim = 100
n = 1000
sphereData = pracma::rands( n, sphereDim, 1. )
mysig = 0.1
spherEmbed = matrix( rnorm( n * embeddDim, 0, mysig ), nrow = n, ncol = embeddDim )
spherEmbed[ , 1:ncol( sphereData ) ] = spherEmbed[ , 1:ncol( sphereData ) ] + sphereData
mysvd = svd( ( cov( spherEmbed ) ) )
plot( ts( mysvd$d[1:22] ) )

Now, we will compute singular values for the covariance matrix, for each radius r and each point on the data set.

calcRowMatDist <- function( xmat, xrow )
  {
  locmag <- function( x, xrow ) sqrt( sum( ( x - xrow )^2 ) )
  apply( xmat, FUN=locmag, MARGIN=1, xrow=xrow )
  }
distmat = ( as.matrix( dist( spherEmbed ) ) )
hist( distmat )

rm( fullCov )

## Warning in rm(fullCov): object 'fullCov' not found

rm( distmat )
nev = sphereDim+10
myp = seq( 0.0, 0.98, 0.02 )
myxaxis = seq( 1.0, 2.2, 0.05 )
mresponse = matrix( ncol = nev, nrow = length( myxaxis ) )
ct = 1
myrct = 0
for ( myr in myxaxis )
  {
  locn = 100
  locsam = sample( 1:n , locn )
  myevs = matrix( nrow=locn, ncol=nev )
  for ( i in 1:locn )
    {
    sel = calcRowMatDist( spherEmbed, spherEmbed[ locsam[i], ] ) < myr
    if ( sum( sel ) >  1 ) {
      lcov = cov( spherEmbed[sel,] )
      temp = svd( lcov )$d[ 1:nev ] # * embeddDim / sum(sel)
#      temp = svd( fullCov[ sel , sel ] )$d[ 1:nev ] * embeddDim / sum(sel)
      } else temp = rep( 0, nev )
    myevs[ i, 1:nev ] = temp
    if ( i == locn ) {
      print( paste( sum(sel), myr ) )
      print(colMeans( myevs, na.rm=T ) )
#      plot( ts( colMeans( myevs, na.rm=T ) ) )
      mresponse[ ct, ] = colMeans( myevs, na.rm=T ) 
      ct = ct + 1
      }
    }
  }

## [1] "1 1"
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1] "1 1.05"
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1] "1 1.1"
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1] "1 1.15"
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1] "1 1.2"
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1] "2 1.25"
##  [1] 1.528638e-02 4.794183e-18 1.526159e-18 1.526159e-18 1.526159e-18
##  [6] 1.526159e-18 1.526159e-18 1.526159e-18 1.526159e-18 1.526159e-18
## [11] 1.526159e-18 1.526159e-18 1.526159e-18 1.526159e-18 1.526159e-18
## [16] 1.526159e-18 1.526159e-18 1.526159e-18 1.526159e-18
## [1] "2 1.3"
##  [1] 5.562819e-02 1.698846e-17 5.493187e-18 5.493187e-18 5.493187e-18
##  [6] 5.493187e-18 5.493187e-18 5.493187e-18 5.493187e-18 5.493187e-18
## [11] 5.493187e-18 5.493187e-18 5.493187e-18 5.493187e-18 5.493187e-18
## [16] 5.493187e-18 5.493187e-18 5.493187e-18 5.493187e-18
## [1] "1 1.35"
##  [1] 2.050274e-01 2.210403e-02 1.918121e-03 1.573832e-03 2.087028e-17
##  [6] 2.045627e-17 2.034627e-17 2.027492e-17 2.009659e-17 2.001469e-17
## [11] 2.001469e-17 2.001469e-17 2.001469e-17 2.001469e-17 2.001469e-17
## [16] 2.001469e-17 2.001469e-17 2.001469e-17 2.001469e-17
## [1] "1 1.4"
##  [1] 3.745064e-01 8.851246e-02 2.756675e-02 1.421336e-02 1.356929e-03
##  [6] 9.698201e-04 7.703651e-04 3.825764e-17 3.691651e-17 3.631206e-17
## [11] 3.601766e-17 3.597265e-17 3.591724e-17 3.580481e-17 3.580374e-17
## [16] 3.574541e-17 3.574541e-17 3.574541e-17 3.574541e-17
## [1] "9 1.45"
##  [1] 4.724290e-01 1.890931e-01 1.055823e-01 5.040993e-02 3.002967e-02
##  [6] 2.092837e-02 1.322148e-02 9.670249e-03 4.742231e-03 3.067700e-03
## [11] 1.939341e-03 4.697706e-04 4.445214e-04 4.753096e-17 4.647612e-17
## [16] 4.583910e-17 4.514188e-17 4.483863e-17 4.449500e-17
## [1] "1 1.5"
##  [1] 0.364134477 0.229484997 0.153101849 0.091682770 0.069707250
##  [6] 0.045743511 0.029918077 0.025469828 0.017512177 0.014043624
## [11] 0.011410220 0.009269545 0.007138760 0.004269362 0.003488454
## [16] 0.002345523 0.002072843 0.001867505 0.001385583
## [1] "12 1.55"
##  [1] 0.261975891 0.204362231 0.162078099 0.127340316 0.109333088
##  [6] 0.084856943 0.065444951 0.054928058 0.043714268 0.035892513
## [11] 0.027477241 0.022489650 0.018274645 0.015095539 0.010917452
## [16] 0.009545192 0.007728845 0.006886075 0.006185826
## [1] "14 1.6"
##  [1] 0.21540697 0.15774438 0.13549265 0.11881491 0.10266237 0.08565662
##  [7] 0.07202900 0.06196917 0.05458017 0.04715279 0.03985059 0.03424481
## [13] 0.02966863 0.02499056 0.02257853 0.01997947 0.01768780 0.01530605
## [19] 0.01371920
## [1] "64 1.65"
##  [1] 0.16453030 0.14048705 0.12376298 0.10843684 0.09648959 0.08222109
##  [7] 0.07073618 0.06280225 0.05439772 0.04568956 0.04106271 0.03681767
## [13] 0.03380111 0.03094558 0.02860257 0.02709419 0.02522771 0.02375017
## [19] 0.02234008
## [1] "36 1.7"
##  [1] 0.14527844 0.12636749 0.11420735 0.10240389 0.09192355 0.08267461
##  [7] 0.07417941 0.06552407 0.05705438 0.04273012 0.03722679 0.03464084
## [13] 0.03271760 0.03076666 0.02903264 0.02736794 0.02575789 0.02449293
## [19] 0.02336848
## [1] "141 1.75"
##  [1] 0.13221838 0.11769696 0.10871860 0.09977901 0.09186515 0.08456593
##  [7] 0.07712654 0.06975966 0.06118922 0.03768633 0.03007532 0.02826100
## [13] 0.02698210 0.02578018 0.02471536 0.02378818 0.02298182 0.02207299
## [19] 0.02139502
## [1] "126 1.8"
##  [1] 0.12987293 0.11904267 0.11066247 0.10150246 0.09433771 0.08705079
##  [7] 0.08060338 0.07423972 0.06539871 0.03749419 0.02728223 0.02572526
## [13] 0.02470374 0.02370256 0.02294284 0.02218533 0.02138799 0.02076027
## [19] 0.02011067
## [1] "165 1.85"
##  [1] 0.12931750 0.11945494 0.11158148 0.10413946 0.09809720 0.09203456
##  [7] 0.08598990 0.07975423 0.07183419 0.03896719 0.02411198 0.02298993
## [13] 0.02205376 0.02132407 0.02066716 0.02014814 0.01949714 0.01903279
## [19] 0.01855782
## [1] "354 1.9"
##  [1] 0.12816415 0.12047758 0.11335518 0.10746091 0.10192453 0.09611720
##  [7] 0.08983651 0.08425143 0.07711710 0.04196170 0.02168297 0.02081803
## [13] 0.02014840 0.01951579 0.01902919 0.01852967 0.01806539 0.01766778
## [19] 0.01728235
## [1] "401 1.95"
##  [1] 0.12770725 0.12138487 0.11565146 0.11030803 0.10459710 0.09999598
##  [7] 0.09473775 0.08955104 0.08285761 0.04734901 0.01998923 0.01926091
## [13] 0.01863712 0.01818935 0.01778390 0.01736140 0.01700020 0.01668447
## [19] 0.01634022
## [1] "536 2"
##  [1] 0.12901448 0.12250199 0.11686289 0.11207581 0.10728818 0.10213506
##  [7] 0.09720854 0.09195521 0.08562361 0.05115812 0.01916373 0.01848539
## [13] 0.01794437 0.01752954 0.01712601 0.01676828 0.01646156 0.01616212
## [19] 0.01586135
## [1] "652 2.05"
##  [1] 0.12957680 0.12385293 0.11786813 0.11291310 0.10889810 0.10420514
##  [7] 0.09912665 0.09403029 0.08830523 0.05670690 0.01852160 0.01789620
## [13] 0.01737176 0.01697140 0.01658012 0.01627933 0.01601275 0.01572088
## [19] 0.01544493
## [1] "761 2.1"
##  [1] 0.12882602 0.12404691 0.11937587 0.11451974 0.11088801 0.10675730
##  [7] 0.10148411 0.09659854 0.09191051 0.06452934 0.01768215 0.01712039
## [13] 0.01665906 0.01630742 0.01596282 0.01567665 0.01540047 0.01516710
## [19] 0.01494065
## [1] "607 2.15"
##  [1] 0.12885222 0.12444423 0.12004169 0.11482315 0.11142233 0.10775991
##  [7] 0.10258957 0.09781142 0.09320748 0.07166144 0.01729877 0.01672116
## [13] 0.01626710 0.01593916 0.01567639 0.01539168 0.01511003 0.01488225
## [19] 0.01466770
## [1] "809 2.2"
##  [1] 0.12828735 0.12445691 0.12008245 0.11517340 0.11207297 0.10894044
##  [7] 0.10423065 0.09895131 0.09492587 0.07899514 0.01698971 0.01646930
## [13] 0.01605698 0.01573487 0.01545549 0.01510535 0.01485968 0.01463897
## [19] 0.01443879

colnames( mresponse ) = paste("EV",1:nev,sep='')
rownames( mresponse ) = paste("Scale",1:length(myxaxis),sep='')
mycols = rainbow( nev )
growthRate1 = magic::shift(mresponse[,1],0)-magic::shift(mresponse[,1],1)*0
plot( myxaxis, growthRate1, type='l', col = mycols[1], main='Evals by scale',
      ylim=c(0.00, max( mresponse[,1]) ), xlab='ball-radius', ylab='Expected Eval' )
for ( i in 2:ncol(mresponse) )
  {
  growthRatek = magic::shift(mresponse[,i],0)-magic::shift(mresponse[,i],1)*0
  points( myxaxis, growthRatek, type='l',col=mycols[i])
  }

Now, similar to above but use the number of nearest k neighbors instead of radius.

rm( fullCov )

## Warning in rm(fullCov): object 'fullCov' not found

distmat = ( as.matrix( dist( spherEmbed ) ) )
nev = sphereDim+3
myxaxis = c( 5:10, 20, 50, 100, 150, 200, 250, 300, 400, 500, 750, 1000 )
myxaxis = c( 10:20 )
mresponse = matrix( ncol = nev, nrow = length( myxaxis ) )
ct = 1
for ( myr in myxaxis )
  {
  sphereLoCov = ANTsR::sparseDistanceMatrix( t(spherEmbed), k = myr, kmetric = "euc" )
  locn = 100 # sample this many local points
  locsam = sample( 1:n , locn )
  myevs = matrix( nrow=locn, ncol=nev )
  for ( i in 1:locn )
    {
    sel = sphereLoCov[, locsam[i] ] > 0
    nevloc = min( c(nev-1,myr-1,sum(sel)-1) )
    if ( sum( sel ) >  0 ) 
      {
      temp = svd( sphereLoCov[ sel , sel ], nv=nevloc )$d[ 1:nevloc ] # * embeddDim / myr
      myevs[ i, 1:nevloc ] = temp
      }
    temp[ is.na( temp ) ] = 0
    if ( i == locn  ) {
      print( paste( sum(sel), myr ) )
      print(colMeans( myevs, na.rm=T ) )
#      plot( ts( colMeans( myevs, na.rm=T ) ) )
      mresponse[ ct, ] = colMeans( myevs, na.rm=T ) 
      ct = ct + 1
      }
    }
  }

## [1] "8 10"
##  [1] 3.6907742 3.3159423 2.3756865 2.2320041 1.3007426 1.1456757 0.6267779
##  [8] 0.6191169 0.6204436       NaN       NaN       NaN
## [1] "10 11"
##  [1] 4.4686170 3.9354013 2.9374966 2.6298473 1.8198039 1.6351505 1.1957445
##  [8] 0.9702918 0.7295170 0.7942155       NaN       NaN
## [1] "3 12"
##  [1] 4.7550441 3.9221084 2.9630362 2.7694864 2.1223764 1.7970147 1.4028699
##  [8] 1.1566263 0.8444663 0.7236365 0.6308166       NaN
## [1] "8 13"
##  [1] 5.964923 4.649441 3.602368 3.275894 2.723695 2.303265 1.980844
##  [8] 1.686827 1.368663 1.248766 1.129746      NaN
## [1] "13 14"
##  [1] 6.376888 4.959728 3.989684 3.611148 2.985291 2.712427 2.206157
##  [8] 1.945518 1.579221 1.367464 1.073171      NaN
## [1] "12 15"
##  [1] 6.990924 5.343704 4.392826 3.922840 3.231498 2.939498 2.497388
##  [8] 2.095609 1.890567 1.604730 1.485234      NaN
## [1] "15 16"
##  [1] 7.009279 5.297982 4.321177 4.011044 3.422107 3.028941 2.720684
##  [8] 2.414744 2.155715 1.893946 1.711060      NaN
## [1] "14 17"
##  [1] 8.434046 5.956580 5.130298 4.704173 4.155646 3.742784 3.294730
##  [8] 2.954401 2.591938 2.318575 2.043359      NaN
## [1] "22 18"
##  [1] 8.399584 6.013943 5.152884 4.588591 4.057728 3.745578 3.299630
##  [8] 2.974994 2.650474 2.385397 2.138050      NaN
## [1] "82 19"
##  [1] 8.880792 6.291354 5.402034 4.888827 4.315951 3.916795 3.451055
##  [8] 3.126454 2.874611 2.585676 2.359225      NaN
## [1] "28 20"
##  [1] 9.852394 6.738374 5.975560 5.474991 4.966476 4.479197 4.068273
##  [8] 3.732722 3.337020 2.992732 2.879458      NaN

colnames( mresponse ) = paste("EV",1:nev,sep='')
rownames( mresponse ) = paste("Scale",1:length(myxaxis),sep='')
mycols = rainbow( nev )
growthRate1 = magic::shift(mresponse[,1],0)-magic::shift(mresponse[,1],1)*0
plot( myxaxis, growthRate1, type='l', col = mycols[1], main='Evals by scale',
      ylim=c(0, max( mresponse[,1]) ), xlab='ball-radius', ylab='Expected Eval' )
for ( i in 2:ncol(mresponse) )
  {
  growthRatek = magic::shift(mresponse[,i],0)-magic::shift(mresponse[,i],1)*0
  points( myxaxis, growthRatek, type='l',col=mycols[i])
  }

sphereLoCov = ANTsR::sparseDistanceMatrix( spherEmbed, k = 50, kmetric = "euc" )
mysvdLoCov = irlba::partial_eigen( scale( sphereLoCov ), n = 5 )
plot( ts( mysvdLoCov$values ) )
# mysvdLoCov = irlba::irlba( sphereLoCov, nv = 10 )
# plot( ts( mysvdLoCov$d ) )

Done!

Blind Source Separation Example

Brian B. Avants

July 14, 2016

blindSourceSeparationInANTsR