# ica <master> <inputFile> <outputFile> <k> <c>
#
# performs ICA
#
# k - number of principal components to use before ica
# c - number of independent components to find
#
# example:
#
# pyspark ica.py local data/sigs.txt results 4 4
#

import sys
import os
from numpy import random, sqrt, zeros, real, dot, outer, diag, transpose, shape
from scipy.linalg import sqrtm, inv, orth
from scipy.io import loadmat
from thunder.util.dataio import parse, saveout
from thunder.factorization.util import svd1, svd2, svd3, svd4
from pyspark import SparkContext
from pyspark.accumulators import AccumulatorParam


class MatrixAccumulatorParam(AccumulatorParam):
    def zero(self, value):
        return zeros(shape(value))

    def addInPlace(self, val1, val2):
        val1 += val2
        return val1

argsIn = sys.argv[1:]
if len(argsIn) < 5:
    print >> sys.stderr, "usage: ica <master> <inputFile> <outputFile> <k> <c>"
    exit(-1)

# parse inputs
sc = SparkContext(argsIn[0], "ica")
dataFile = str(argsIn[1])
outputDir = str(argsIn[2]) + "-ica"
k = int(argsIn[3])
c = int(argsIn[4])
if not os.path.exists(outputDir):
    os.makedirs(outputDir)

# load data
lines = sc.textFile(dataFile)
data = parse(lines, "raw", None, [150, 1000]).cache()
n = data.count()

# reduce dimensionality
#comps, latent, scores = svd4(sc, data, k, 0)

# whiten data
#whtMat = real(dot(inv(diag(sqrt(latent))), comps))
#unwhtMat = real(dot(transpose(comps), diag(sqrt(latent))))

whtMat = loadmat(outputDir + "/whtMat.mat")['whtMat']
unwhtMat = loadmat(outputDir + "/unwhtMat.mat")['unwhtMat']
wht = data.map(lambda x: dot(whtMat, x))
#print(wht.first())

# save whitening matrices
#saveout(whtMat, outputDir, "whtMat", "matlab")
#saveout(unwhtMat, outputDir, "unwhtMat", "matlab")

# do multiple independent component extraction
B = orth(random.randn(k, c))
#B = loadmat(outputDir + "/B.mat")['B']
Bold = zeros((k, c))
iterNum = 0
minAbsCos = 0
tol = 0.000001
iterMax = 1000
errVec = zeros(iterMax)


global Bnew

def outerSum(iterator, B):
    yield sum(outer(x, dot(x, B) ** 3) for x in iterator)

def matrixAccum(x):
    global Bnew
    Bnew += x

def outerSumAccum(x):
    global Bnew
    Bnew += outer(x, x)

def outerSumAccumOther(x, B):
    global Bnew
    Bnew += outer(x, dot(x, B) ** 3)

while (iterNum < iterMax) & ((1 - minAbsCos) > tol):
    iterNum += 1
    # update rule for pow3 nonlinearity (TODO: add other nonlins)
    #Bnew = sc.accumulator(zeros((k, c)), MatrixAccumulatorParam())
    #wht.foreach(lambda x: outerSumAccumOther(x, B))
    #wht.mapPartitions(lambda x: outerSum(x, B)).foreach(matrixAccum)
    #wht.foreach(lambda x: outerSumAccumOther(x, B))
    #B = Bnew.value / n - 3 * B
    B = wht.map(lambda x: outer(x, dot(x, B) ** 3)).sum() / n - 3 * B
    # orthognalize
    B = dot(B, real(sqrtm(inv(dot(transpose(B), B)))))
    # evaluate error
    minAbsCos = min(abs(diag(dot(transpose(B), Bold))))
    # store results
    Bold = B
    saveout(B, outputDir, "B", "matlab")
    errVec[iterNum-1] = (1 - minAbsCos)

# get unmixing matrix
W = dot(transpose(B), whtMat)

# get components
sigs = data.map(lambda x: dot(W, x))

# save output files
saveout(W, outputDir, "W", "matlab")
saveout(sigs, outputDir, "sigs", "matlab", c)

@@ -547,249 +547,8 @@
 ext%0A
-from pyspark.accumulators import AccumulatorParam%0A%0A%0Aclass MatrixAccumulatorParam(AccumulatorParam):%0A    def zero(self, value):%0A        return zeros(shape(value))%0A%0A    def addInPlace(self, val1, val2):%0A        val1 += val2%0A        return val1
 %0A%0Aar
@@ -1556,16 +1556,17 @@
 raction%0A
+#
 B = orth
@@ -1586,17 +1586,16 @@
 (k, c))%0A
-#
 B = load
@@ -1732,716 +1732,143 @@
 )%0A%0A%0A
-global Bnew%0A%0Adef outerSum(iterator, B):%0A    yield sum(outer(x, dot(x, B) ** 3) for x in iterator)%0A%0Adef matrixAccum(x):%0A    global Bnew%0A    Bnew += x%0A%0Adef outerSumAccum(x):%0A    global Bnew%0A    Bnew += outer(x, x)%0A%0Adef outerSumAccumOther(x, B):%0A    global Bnew%0A    Bnew += outer(x, dot(x, B) ** 3)%0A%0Awhile (iterNum %3C iterMax) & ((1 - minAbsCos) %3E tol):%0A    iterNum += 1%0A    # update rule for pow3 nonlinearity (TODO: add other nonlins)%0A    #Bnew = sc.accumulator(zeros((k, c)), MatrixAccumulatorParam())%0A    #wht.foreach(lambda x: outerSumAccumOther(x, B))%0A    #wht.mapPartitions(lambda x: outerSum(x, B)).foreach(matrixAccum)%0A    #wht.foreach(lambda x: outerSumAccumOther(x, B))%0A    #B = Bnew.value / n - 3 * B
+while (iterNum %3C iterMax) & ((1 - minAbsCos) %3E tol):%0A    iterNum += 1%0A    # update rule for pow3 nonlinearity (TODO: add other nonlins)
 %0A   

Use previous B