hnfa.py

Extension packages to Bayes Blocks library, reported in "Nonlinear independent factor analysis by hi

            vA_out = self.net.GetGaussianNode(Label("vA_out",
                                                    1, self.nexthidden))
            prod = []
            for i in range(len(weights[j])):
                if weights[j][i] is None:
                    continue
                A = self.fact.GetGaussian(
		    Label("A", 2, 1, i, self.nexthidden),
		    self.fact.GetConst0(), vA_in)
                self.fact.EvidenceNode(
                    A, mean=weights[j][i],
                    var=0.01, decay=40.0)
                prod.append(self.fact.GetProdV(
			Label("prA", 2, 1, i, self.nexthidden),
			A, self.net.GetNode(Label("s", 2, i))))
                A.SetPersist(5)
            ms = (self.fact.GetGaussian(Label("ms", 1, self.nexthidden),
                               self.fact.GetConst0(),
                               self.fact.GetConst0()))
            ms.SetPersist(5)
            self.fact.EvidenceNode(
                ms, mean=bias[j],
                var=0.01, decay=40.0)
            sum = self.fact.BuildSum2VTree(prod + [ms], ("sums",
                                                        1, self.nexthidden))
            vs = self.net.GetGaussianNode(Label("vs", 1, self.nexthidden))
            self.fact.EvidenceNode(vs, mean=9, var=1, decay=vsdecay)
            s = self.fact.GetGaussNonlinV(Label("s", 1, self.nexthidden),
					  sum, vs)
            s.Update()
	    s.SetPersist(5)
            for i in childind[j]:
                vA = self.fact.GetSum2(Label("vA", 1, 0, self.nexthidden, i),
                          vA_out, self.net.GetGaussianNode(Label("vA_in", 0, i)))
                A = self.fact.GetGaussian(Label("A", 1, 0, self.nexthidden, i),
					  self.fact.GetConst0(), vA)
                prA = self.fact.GetProdV(
		    Label("prA", 1, 0, self.nexthidden, i),
		    A, s)
                sumcache[i].append(prA)
                Alist.append(A)
            self.nexthidden += 1
        for i in range(self.dim[0]):
            if len(sumcache[i]) == 0:
                continue
            sum = self.fact.BuildSum2VTree(sumcache[i], ("sums", 0, i))
            Helpers.AddSum2V(self.net,
                             self.net.GetNode(Label("s", 0, i)),
                             0, sum, Label("sums", 0, i))
        map(Net.Variable.Update,Alist)
        self.net.SortNodes()

    def AddHiddenBest(self, num=5, numtest=1000,
                      scalemean = 0.5, scalestd = 0.5, vsdecay = 1000):
        if num == numtest:
            return self.AddHidden(num=num,scalemean=scalemean,
                                  scalestd=scalestd,vsdecay=vsdecay)
        (weights, bias, costchange)=self.HiddenCost(
            numtest=numtest, scalemean=scalemean, scalestd=scalestd)
        cs = map(Numeric.sum, costchange)
        s=list(Numeric.argsort(cs))
        ind = s[-5:]
        ind.reverse()
        self.AddHidden2(weights=Numeric.take(weights,ind),
                        bias=Numeric.take(bias,ind),
                        vsdecay=vsdecay)

    def HiddenCost(self, numtest=1, scalemean = 0.5, scalestd = 0.5):
        sources = Numeric.array(Helpers.GetMeanV(
            map(lambda i:self.net.GetVariable(Label('s', 2, i)),
                range(self.dim[2]))),'d')
        gradmean = Numeric.zeros((self.dim[0], self.net.Time()),'d')
        gradvar = Numeric.zeros((self.dim[0], self.net.Time()),'d')
        exp = []
        for i in range(self.dim[0]):
            s0 = self.net.GetVariable(Label('s', 0, i))
            dsset = Net.DVSet()
            s0.GradRealV(dsset, s0.GetParent(0))
            gradmean[i] = dsset.mean
            gradvar[i] = dsset.var
            exp.append(Helpers.GetExp(s0.GetParent(1)))

        weights=[]
        bias=[]
        costchange=[]
        for j in range(numtest):
            scale = math.exp(self.random.gauss(scalemean, scalestd))
            weights.append([])
            for i in range(self.dim[2]):
                weights[j].append(self.random.gauss(
                    0,scale/math.sqrt(self.dim[2])))
            bias.append(self.random.gauss(0, scale))
            try:
                mean1=Numeric.dot(Numeric.array([weights[j]]), sources)+bias[j]
                #This clip should take care of overflows
                mean1=Numeric.clip(mean1,-18,18)
                mean2=Numeric.exp(-(mean1**2))
                mean2=mean2-MLab.mean(mean2,1)
                #mean2=Numeric.exp(-(mean1**2)/(1+2*mvs1))/math.sqrt(1+2*mvs1)
                #Want to minimise cost2 = sum(a*(x-b*m))^2
                #where a=a[i], x=x[i], m=m[i], b is the parameter
                #gradmean=a*x gradvar=a/2 mean2=m
                #differentiate-> sum(-m*a*(x-b*m))==0
                #-> b=sum(m*a*x)/sum(m^2*a) = c/d
                #we wan the difference sum(a*x^2-a*(x-b*m)^2)=
                #sum(2*a*x*b*m-a*b^2*m^2)=
                #2*b*sum(a*x*m)-b^2*sum(a*m^2)=2*b*c-b^2*d=2*c^2/d-c^2/d=c^2/d
                c=Numeric.dot(mean2,Numeric.transpose(gradmean))
                _2d=Numeric.dot(mean2**2, Numeric.transpose(gradvar))
                costchange.append(((c**2)/_2d)[0])
            except OverflowError:
                costchange.append(None)
        return weights, bias, costchange

    def PruneHidden(self, hidden, cdiff = None):
        node = self.net.GetNode(hidden)
        if cdiff is None:
            cdiff = Helpers.CostDifferenceV(node)
        self.learner.HistoryAdd("PruneHidden", (hidden, cdiff))
        node.Die()
        self.net.CleanUp()
        self.net.SortNodes()

    def TryPruneHidden(self, cdiff = 200, num=1, lazy=0):
        self.learner.HistoryAdd("TryPruneHidden", (cdiff, num))
        nodes = filter(lambda x:Unlabel(x.GetLabel())[1][1] <
                       self.nexthidden-lazy,
                       self.net.GetVariables('s\(1,'))
        hiddens = map(lambda x:(Helpers.CostDifferenceV(x), x.GetLabel()),
                      nodes)
        hiddens.sort()
        for i in range(num):
            if len(hiddens) > i and hiddens[i][0] < cdiff:
                self.PruneHidden(hiddens[i][1], hiddens[i][0])
            else:
                return i
        return num

    def AddWeight(self):
        raise RuntimeError, "AddWeight() not ready for use"
        hiddenlist=self.net.GetVariables('s\(1,')
        hidden=self.random.choise(slist)
        hiddennum = Unlabel(s.GetLabel())[1][1]
        parentind=self.random.randrange(self.dim[2])
        childind=self.random.randrange(self.dim[0])
        vA_in = self.net.GetGaussianNode(Label("vA_in", 1, hiddennum))
        vA_out = self.net.GetGaussianNode(Label("vA_out", 1, hiddennum))

    def AddAllWeights(self, prob = 1.0):
        adden = []
        for i in range(self.dim[0]):
            for j in range(self.dim[2]):
                if self.random.random() < prob:
                    if self.AddOneWeight((2, 0, j ,i)):
                        adden.append(Label('A', 2, 0, j ,i))
        hiddenlist=self.net.GetVariables('s\(1,')
        for h in hiddenlist:
            k = Unlabel(h.GetLabel())[1][1]
            for j in range(self.dim[2]):
                if self.random.random() < prob:
                    if self.AddOneWeight((2, 1, j ,k)):
                        adden.append(Label('A', 2, 1, j ,k))
            for i in range(self.dim[0]):
                if self.random.random() < prob:
                    if self.AddOneWeight((1, 0, k ,i)):
                        adden.append(Label('A', 1, 0, k ,i))
        self.net.SortNodes()
        self.learner.HistoryAdd("AddAllWeights", (prob, tuple(adden)))
            

    def AddOneWeight(self, index, lazy = 0):
        if self.net.GetNode(Label("A", index)) is None:
            if index[0] == 1:
                vA_out = self.net.GetGaussianNode(Label("vA_out", 1, index[2]))
            elif index[0] == 2:
                vA_out = None
            else:
                raise IndexError, "index[0] out of range"
            if index[1] not in (0, 1):
                raise IndexError, "index[1] out of range"
            vA_in = self.net.GetGaussianNode(Label("vA_in", index[1], index[3]))
            if vA_out is None:
                vA = vA_in
            else:
                vA = self.fact.GetSum2(Label("vA", index), vA_out, vA_in)
            A = self.fact.GetGaussian(Label("A", index),
				      self.fact.GetConst0(), vA)
            A.SetPersist(5)
            prod = self.fact.GetProdV(Label("prA", index), A,
                             self.net.GetNode(Label("s", index[0], index[2])))
            if lazy:
                return prod
            Helpers.AddSum2V(self.net,
                             self.net.GetNode(Label("s",
                                                    index[1], index[3])),
                             0, prod,
                             Label("sums", index[1], index[3]))
            A.Update()
            self.net.SortNodes()
            return 1

    def TryPruning(self):
        return self.learner.TryPruning()

    def HookeJeeves(self):
        return self.learner.HookeJeeves()
        
    def Iteration(self):
        return self.learner.Iteration()

    def LearnNet(self, *args, **kws):
        return self.learner.LearnNet(*args, **kws)