# HG changeset patch
# User Joseph Turian <turian@gmail.com>
# Date 1236582680 14400
# Node ID 12b1b09ffd2b32349b0c507331b736d00c0453de
# Parent  85436cda77baa08525d7a1195506c14e717b14f0
Added preliminary code for computing negative log Poisson cost

diff -r 85436cda77ba -r 12b1b09ffd2b pylearn/algorithms/sandbox/cost.py
--- a/pylearn/algorithms/sandbox/cost.py	Mon Mar 09 00:25:46 2009 -0400
+++ b/pylearn/algorithms/sandbox/cost.py	Mon Mar 09 03:11:20 2009 -0400
@@ -45,9 +45,12 @@
 def nlpoisson(target, output, beta_scale=1, axis=0):
     """
     The negative log Poisson regression probability.
-    From Marc'Aurelio and Szummer (2008).
+    From Ranzato and Szummer (2008).
 
     Output should be of the form Weight*code+bias, i.e. unsquashed.
+    NB this is different than the formulation in Salakhutdinov and Hinton
+    (2007), in which the output is softmax'ed and multiplied by the
+    input document length.
 
     There is a beta term that is proportional to document length. We
     are not sure what beta scale is used by the authors. We use 1 as
@@ -55,7 +58,37 @@
 
     Axis is the axis along which we sum the target values, to obtain
     the document length.
-    @bug: This axis may be wrong!!
     """
-    beta = beta_scale * T.sum(target, axis=axis)
-    return beta * T.exp(output) - T.dot(target, output) + logfactorial(target)
+#    from theano.printing import Print
+    doclen = tensor.sum(target, axis=axis)
+    beta = beta_scale * doclen
+    return tensor.sum(beta * tensor.exp(output) - target*output + logfactorial(target), axis=axis)
+
+
+#import numpy
+#def nlpoisson_nontheano(target, output, beta_scale=1, axis=0):
+#    doclen = numpy.sum(target, axis=axis)
+#    print "doclen", doclen
+#    beta = beta_scale * doclen
+#    print "beta", beta
+#    print "exp", numpy.exp(output)
+#    print "beta * exp", beta * numpy.exp(output)
+#    print "x * y", target * output
+#
+#    import theano.tensor as TT
+#    x = TT.as_tensor(target)
+#    o = logfactorial(x)
+#    f = T.function([],o)
+#    logf = f()
+#    print "log factorial(x)", logf
+#    print "beta * exp - dot + log factorial", beta * numpy.exp(output) - target*output + f()
+#    print "total loss", numpy.sum(beta * numpy.exp(output) - target*output + f(), axis=axis)
+#
+##    return beta * numpy.exp(output) - numpy.dot(target, output)
+##            #+ logfactorial(target)
+#
+#import numpy
+#target = numpy.array([0, 0, 1, 1, 2, 2, 100, 100])
+##output = numpy.array([0., 0.5, 1., 0.5, 2., 0.5, 100., 0.5])
+#output = numpy.array([0., 1, 1., 0, 1, 0, 5, 1])
+#nlpoisson_nontheano(target, output)
diff -r 85436cda77ba -r 12b1b09ffd2b pylearn/algorithms/sandbox/test_cost.py
--- a/pylearn/algorithms/sandbox/test_cost.py	Mon Mar 09 00:25:46 2009 -0400
+++ b/pylearn/algorithms/sandbox/test_cost.py	Mon Mar 09 03:11:20 2009 -0400
@@ -22,5 +22,20 @@
 #        print repr(f())
         self.failUnless(numpy.all(f() == numpy.asarray([0., 0., 1.38629436, 3.29583687, 5.54517744, 8.04718956, 10.75055682, 13.62137104, 16.63553233, 19.7750212])))
 
+class T_nlpoisson(unittest.TestCase):
+    def test(self):
+        target = TT.as_tensor([0, 0, 1, 1, 2, 2, 100, 100])
+        output = TT.as_tensor([0., 1, 1., 0, 1, 0, 5, 1])
+        o = cost.nlpoisson(target, output)
+        f = T.function([],o)
+        self.failUnless(f() - 33751.7816277 < 1e-5)
+
+#    def test_gradient(self):
+#        target = TT.as_tensor([0, 0, 1, 1, 2, 2, 100, 100])
+#        output = TT.as_tensor([0., 1, 1., 0, 1, 0, 5, 1])
+#        o = cost.nlpoisson(target, output)
+#        f = T.function([],o)
+#        self.failUnless(f() - 33751.7816277 < 1e-5)
+
 if __name__ == '__main__':
     unittest.main()