m2 = nf.GetGaussian("m2", c0, c_5)        v2 = nf.GetGaussian("v2", c0, c_5)        nf.EvidenceNode(m1, mean=-1.0)        nf.EvidenceNode(m2, mean=1.0)        vx = nf.GetConstant("vx", 9.21)        d = []        s = []        x = []        for i in range(dim):            d.append(nf.GetDiscreteDirichlet(Label("d", i), c))            s.append(nf.GetMoG(Label("s", i), d[i]))            s[i].AddComponent(m1, v1)            s[i].AddComponent(m2, v2)            x.append(nf.GetGaussian(Label("x", i), s[i], vx))            x[i].Clamp(self.data[i])        return net    def costdump(self, vector, scalar):        variables = ["c", "m1", "v1", "m2", "v2"]        for l in variables:            n1 = vector.GetVariable(l)            n2 = scalar.GetVariable(l)            print "%s : %f : %f" % (l, n1.Cost(), n2.Cost())        variables = ["d", "s", "x"]        for l in variables:            c1 = vector.GetVariable(l).Cost()            c2 = MLab.sum([n.Cost() for n in scalar.GetVariableArray(l)])            print "%s : %f : %f" % (l, c1, c2)    def comparenets(self, vector, scalar):        dim = self.data.shape[0]        #self.costdump(vector, scalar)        # this alone is quite reliable indicator        costdiff = abs(vector.Cost() - scalar.Cost())        self.assert_(costdiff < dim * 1e-8, "costdiff = %f" % costdiff)        s1 = GetVarV(vector.GetVariable("s"))        s2 = GetVar(scalar.GetVariableArray("s"))####        import plot##        pl = plot.Plotter()##        pl.multiplot(Numeric.log10([s1,s2]))        self.assert_(MLab.max(Numeric.absolute(s1 - s2)) < 1e-8)        for l in ["m1", "v1", "m2", "v2"]:            n1 = GetMean(vector.GetVariable(l))            n2 = GetMean(scalar.GetVariable(l))            #print "%s : %f : %f" % (l, n1, n2)            self.assert_(abs(n1 - n2) < 1e-8)        d1 = GetDiscreteV(vector.GetVariable("d"))        d2 = Numeric.transpose(GetDiscrete(scalar.GetVariableArray("d")))        #print d1        #print d2        e = MLab.max(Numeric.reshape(            Numeric.absolute(d1 - d2), (d1.shape[0]*d1.shape[1],)))        self.assert_(e < 1e-8, "max e = %f" % e)        c1 = vector.GetVariable("c")        c2 = scalar.GetVariable("c")        #print GetMeanV(c1)        self.assert_(            MLab.mean(Numeric.absolute(GetMeanV(c1) - GetMeanV(c2))) < 1e-8)        self.assert_(            MLab.mean(Numeric.absolute(GetVarV(c1) - GetVarV(c2))) < 1e-8)            def test1(self):        """Sameness of scalar and vector implementation of MoGs."""        net1 = self.buildnet1()        net2 = self.buildnet2()                for i in range(50):            net1.UpdateAll()            net2.UpdateAll()            #print "%d : %f : %f" % (i, net1.Cost(), net2.Cost())        self.comparenets(net1, net2)    def test2(self):        """Load/save interface."""        net1 = self.buildnet1()        net2 = self.buildnet2()        for i in range(50):            net1.UpdateAll()            net2.UpdateAll()        self.comparenets(net1, net2)        from PickleHelpers import load, save        save(net2, "/tmp/bblocks_tmpnet")        net2alt = load("/tmp/bblocks_tmpnet")        self.comparenets(net1, net2alt)class GaussRectTestCase(unittest.TestCase):    def setUp(self):        t = Numeric.arange(20) / float(20)        s = Numeric.sin(4*Numeric.pi*t) + 0.5        # N(0, 0.01^2)        noise = Numeric.array(            [-0.00020618, -0.00214195, -0.00860223,  0.0083709 ,             0.00493691, -0.00303564,  0.00467232,  0.00875158,             0.01462178,  0.00804469,  0.00791851,  0.00444479,             -0.00398483, -0.00677404, -0.02249393,  0.00172111,            -0.00420562,  0.00287512,  0.02004725,  0.00479894],            Numeric.Float)        self.data = s * (s > 0) + noise    def buildstatnet1(self):        net = PyNet.PyNet(self.data.shape[0])        nf = PyNet.PyNodeFactory(net)        c0 = nf.GetConstant("c0", 0.0)        c_5 = nf.GetConstant("c_5", -5.0)        ms = nf.GetGaussian("ms", c0, c_5)        vs = nf.GetGaussian("vs", c0, c_5)        nf.EvidenceNode(ms, mean=1.0)        nf.EvidenceNode(vs, mean=0.0)        s = nf.GetGaussRectV("s", ms, vs)        r = nf.GetRectificationV("r", s)        vx = nf.GetConstant("vx", -2*Numeric.log(0.01))        #vx = nf.GetGaussian("vx", c0, c_5)        x = nf.GetGaussianV("x", r, vx)        x.Clamp(Helpers.Array2DV(self.data))        return net    def buildstatnet2(self):        net = PyNet.PyNet(1)        nf = PyNet.PyNodeFactory(net)        c0 = nf.GetConstant("c0", 0.0)        c_5 = nf.GetConstant("c_5", -5.0)        ms = nf.GetGaussian("ms", c0, c_5)        vs = nf.GetGaussian("vs", c0, c_5)        nf.EvidenceNode(ms, mean=1.0)        nf.EvidenceNode(vs, mean=0.0)        vx = nf.GetConstant("vx", -2*Numeric.log(0.01))        for i in range(self.data.shape[0]):            s = nf.GetGaussRect(Label("s", i), ms, vs)            r = nf.GetRectification(Label("r", i), s)            x = nf.GetGaussian(Label("x", i), r, vx)            x.Clamp(self.data[i])        return net    def comparenets(self, vector, scalar):        dim = self.data.shape[0]        costdiff = abs(vector.Cost() - scalar.Cost())        self.assert_(costdiff < dim * 1e-8, "costdiff = %f" % costdiff)        s1 = GetMeanV(vector.GetVariable("s"))        s2 = GetMean(scalar.GetVariableArray("s"))        e = MLab.max(Numeric.absolute(s1 - s2))        self.assert_(e < 1e-8)                r1 = GetMeanV(vector.GetNode("r"))        r2 = GetMean(scalar.GetNodeArray("r"))                e = MLab.max(Numeric.absolute(r1 - r2))        self.assert_(e < 1e-8)        if vector.GetVariable("ms") is not None:            self.assert_(abs(GetMean(vector.GetVariable("ms"))                             - GetMean(scalar.GetVariable("ms"))) < 1e-8)        else:            self.assert_(abs(GetMean(vector.GetVariable("m0"))                             - GetMean(scalar.GetVariable("m0"))) < 1e-8)        self.assert_(abs(GetMean(vector.GetVariable("vs"))                         - GetMean(scalar.GetVariable("vs"))) < 1e-8)            def test1(self):        """Vector/scalar comparison with static model."""        net1 = self.buildstatnet1()        net2 = self.buildstatnet2()                for i in range(50):            net1.UpdateAll()            net2.UpdateAll()                        #print "%d : %f : %f" % (i, net1.Cost(), net2.Cost())        self.comparenets(net1, net2)    def builddynnet1(self):        net = PyNet.PyNet(self.data.shape[0])        nf = PyNet.PyNodeFactory(net)        c0 = nf.GetConstant("c0", 0.0)        c_5 = nf.GetConstant("c_5", -5.0)        m0 = nf.GetGaussian("ms", c0, c_5)        vs = nf.GetGaussian("vs", c0, c_5)        proxy = nf.GetProxy("proxy", "s")        delay = nf.GetDelayV("delay", m0, proxy)        s = nf.GetGaussRectV("s", delay, vs)        r = nf.GetRectificationV("r", s)        vx = nf.GetConstant("vx", -2*Numeric.log(0.01))        #vx = nf.GetGaussian("vx", c0, c_5)        x = nf.GetGaussianV("x", r, vx)        x.Clamp(Helpers.Array2DV(self.data))        net.ConnectProxies()        return net    def builddynnet2(self):        net = PyNet.PyNet(1)        nf = PyNet.PyNodeFactory(net)        c0 = nf.GetConstant("c0", 0.0)        c_5 = nf.GetConstant("c_5", -5.0)        m0 = nf.GetGaussian("ms", c0, c_5)        vs = nf.GetGaussian("vs", c0, c_5)        vx = nf.GetConstant("vx", -2*Numeric.log(0.01))        prev = m0        for i in range(self.data.shape[0]):            s = nf.GetGaussRect(Label("s", i), prev, vs)            prev = s            r = nf.GetRectification(Label("r", i), s)            x = nf.GetGaussian(Label("x", i), r, vx)            x.Clamp(self.data[i])        return net    def test2(self):        """Vector/scalar comparison with dynamic model."""        net1 = self.builddynnet1()        net2 = self.builddynnet2()                for i in range(100):            net1.UpdateAll()            net2.UpdateAll()                        #print "%d : %f : %f" % (i, net1.Cost(), net2.Cost())        self.comparenets(net1, net2)    def test3(self):        """Load/save vector."""        from PickleHelpers import load, save                ref = self.buildstatnet2()        net = self.buildstatnet1()        for i in range(100):            ref.UpdateAll()            net.UpdateAll()        self.comparenets(net, ref)        save(net, "/tmp/bblocks_tmpnet")        alt = load("/tmp/bblocks_tmpnet")        self.comparenets(alt, ref)        ref = self.builddynnet2()        net = self.builddynnet1()        for i in range(100):            ref.UpdateAll()            net.UpdateAll()        self.comparenets(net, ref)        save(net, "/tmp/bblocks_tmpnet")        alt = load("/tmp/bblocks_tmpnet")        self.comparenets(alt, ref)            def test4(self):        """Load/save scalar."""        from PickleHelpers import load, save        ref = self.buildstatnet1()        net = self.buildstatnet2()        for i in range(100):            ref.UpdateAll()            net.UpdateAll()        self.comparenets(ref, net)        save(net, "/tmp/bblocks_tmpnet")        alt = load("/tmp/bblocks_tmpnet")        self.comparenets(ref, alt)        ref = self.builddynnet1()        net = self.builddynnet2()        for i in range(100):            ref.UpdateAll()            net.UpdateAll()        self.comparenets(ref, net)        save(net, "/tmp/bblocks_tmpnet")        alt = load("/tmp/bblocks_tmpnet")        self.comparenets(ref, alt)if __name__ == '__main__':    cases = [s for s in dir() if s.endswith("TestCase")]    suites = []    if len(sys.argv) > 1:        for i in range(1, len(sys.argv)):            if sys.argv[i] in cases:                suites.append(unittest.makeSuite(eval(sys.argv[i]), "test"))            else:                print sys.argv[i] + "?"    else:        for case in cases:            suites.append(unittest.makeSuite(eval(case), "test"))    bigsuite = unittest.TestSuite(suites)    runner = unittest.TextTestRunner()    runner.run(bigsuite)##    import pdb##    pdb.run("runner.run(bigsuite)")##