提交 93b4e940 authored 作者: Olivier Breuleux's avatar Olivier Breuleux

cleaned up the files

上级 c9221f08
import unittest
from wrappers import *
class _testCase_input(unittest.TestCase):
def setUp(self):
literal.hdb = {}
literal.udb = {}
def test_input_int(self):
w = input(3)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.int_dtype)
self.failUnless(w.data == 3)
def test_input_float(self):
w = input(3.0)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.float_dtype)
self.failUnless(w.data == 3.0)
class _testCase_wrap(unittest.TestCase):
def setUp(self):
literal.hdb = {}
literal.udb = {}
def test_wrap_int(self):
w = wrap(3)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.int_dtype)
self.failUnless(w.data == 3)
def test_wrap_float(self):
w = wrap(3.0)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.float_dtype)
self.failUnless(w.data == 3.0)
class _testCase_literal(unittest.TestCase):
def setUp(self):
literal.hdb = {}
literal.udb = {}
def test_int(self):
w = literal(3)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.int_dtype)
self.failUnless(w.data == 3)
u = literal(1+2)
self.failUnless(u is w)
def test_float(self):
w = literal(3.0)
self.failUnless(isinstance(w, input.NN))
self.failUnless(str(w.data.dtype) == input.float_dtype)
self.failUnless(w.data == 3.0)
u = literal(1.0+2.0)
self.failUnless(u is w)
def test_mixed(self):
f = literal(2.0)
i = literal(2)
self.failUnless(i is not f)
if __name__ == '__main__':
unittest.main()
import unittest
import constructor_fodder as cf
class Allocator:
def __init__(self, cls, ctor):
self.cls = cls
self.ctor = ctor
def __call__(self, *args, **kwargs):
rval = self.cls.__new__(self.cls, *args, **kwargs)
rval.__init__(*args, **kwargs)
return rval
class ModeOpAllocator:
def __init__(self, cls, ctor):
self.cls = cls
self.ctor = ctor
def __call__(self, *args, **kwargs):
op = self.cls.__new__(self.cls, *args, **kwargs)
op.__init__(*args, **kwargs)
mode = self.ctor.mode()
if mode == 'eval':
op.perform()
if op.nout == 1:
return op.out.data
else:
return [output.data for output in op.outputs]
elif mode == 'build_eval':
op.perform()
if op.nout == 1:
return op.out
else:
return op.outputs
class Constructor:
def __init__(self):
pass
def add_module(self, module, module_name, accept=lambda x:issubclass(x, cf.base)):
dct = {}
for symbol in dir(module):
if symbol[:2] == '__': continue
obj = getattr(module, symbol)
if accept(obj): dct[symbol] = Allocator(obj)
class Dummy:pass
self.__dict__[module_name] = Dummy()
self.__dict__[module_name].__dict__.update(dct)
def add_from_module(self, module, accept=lambda x:issubclass(x, cf.base)):
for symbol in dir(module):
if symbol[:2] == '__': continue
obj = getattr(module, symbol)
#print 'considering', symbol, obj
if accept(obj): self.__dict__[symbol] = Allocator(obj)
def add_globals_from_module(self, module, accept=lambda x:issubclass(x, cf.base)):
for symbol in dir(module):
if symbol[:2] == '__': continue
obj = getattr(module, symbol)
#print 'considering', symbol, obj
if accept(obj):
if hasattr(globals(), symbol):
print 'Warning, overwriting global variable: %s' % symbol
globals()[symbol] = Allocator(obj)
if __name__=='__main__':
c = Constructor()
c.add_module(cf,'cf')
aa,bb = c.cf.A(), c.cf.B()
print aa,bb
c.add_from_module(cf)
a,b = c.A(), c.B()
print a,b
c.add_globals_from_module(cf)
d,e = A(), B()
print d,e
class base(object): pass
class A(base): pass
class B(base): pass
class C(base): pass
import os # for building the location of the .omega/omega_compiled cache directory
import sys # for adding the inline code cache to the include path
import platform #
import unittest
import weakref
import inspect
import md5
import copy
import numpy
from scipy import weave
import gof
from gof import current_mode, set_mode, build_mode, eval_mode, build_eval_mode
from gof import pop_mode, is_result, ResultBase
import type_spec
import cutils
import blas
import compile
# __all__ = ['set_mode', 'get_mode', 'NumpyR', 'NumpyOp']
def build(f, *args, **kwargs):
build_mode()
r = f(*args, **kwargs)
pop_mode()
return r
def as_string(*rs):
s = gof.graph.as_string(gof.graph.inputs(rs), rs)
if len(rs) == 1:
return s[1:-1]
else:
return s
def print_graph(*rs):
print as_string(*rs)
#useful mostly for unit tests
def _approx_eq(a,b,eps=1.0e-9):
a = numpy.asarray(a)
b = numpy.asarray(b)
if a.shape != b.shape:
return False
return numpy.max(numpy.abs(a-b)) < eps
# This function is only executed the first time it is called, subsequent calls
# return immediately from a cache of the first return value
@blas._constant # TODO: move this decorator to a utility file
def _compile_dir():
"""Return the directory in which scipy.weave should store code objects.
If the environment variable OMEGA_COMPILEDIR is set, its value is returned.
If not, a directory of the form $HOME/.omega/compiledir_<platform Id>.
As a test, this function touches the file __init__.py in the returned
directory, and raises OSError if there's a problem.
A directory coming from OMEGA_COMPILEDIR is not created automatically, but
a directory in $HOME/.omega is created automatically.
This directory is appended to the sys.path search path before being
returned, if the touch was successful.
"""
if os.getenv('OMEGA_COMPILEDIR'):
cachedir = os.getenv('OMEGA_COMPILEDIR')
else:
# use (and possibly create) a default code cache location
platform_id = platform.platform() + '-' + platform.processor()
import re
platform_id = re.sub("[\(\)\s]+", "_", platform_id)
cachedir = os.path.join(os.getenv('HOME'), '.omega', 'compiledir_'+platform_id)
if not os.access(cachedir, os.R_OK | os.W_OK):
#this may raise a number of problems, I think all of which are serious.
os.makedirs(cachedir, 7<<6)
cachedir_init = cachedir+'/__init__.py'
touch = os.system('touch '+cachedir_init)
if touch:
raise OSError('touch %s returned %i' % (cachedir_init, touch))
if cachedir not in sys.path:
sys.path.append(cachedir)
return cachedir
class Numpy2(ResultBase):
"""Result storing a numpy ndarray"""
__slots__ = ['_dtype', '_shape', '_order']
class ShapeUnknown: pass # TODO: use this as the shape of uncomputed ndarrays of unknown shape
class StateError(Exception): pass
def __init__(self, role=None, data=None, constant=False):
self._order = 'C'
if isinstance(data, (tuple, list)): # unallocated setup
shape, dtype = data
ResultBase.__init__(self, role, data=None, constant=constant)
self._shape = shape
self._dtype = dtype
else: # allocated setup
ResultBase.__init__(self, role, data, constant)
################################
# ResultBase
#
def data_filter(self, data):
return numpy.asarray(data)
################################
# Numpy2 specific functionality
#
__array__ = property(lambda self: self.data.__array__)
__array_struct__ = property(lambda self: self.data.__array_struct__)
def data_alloc(self):
return numpy.ndarray(shape=self.shape, dtype=self.dtype, order=self._order)
# self._dtype is used when self.data hasn't been set yet
def __dtype_get(self):
if self.data is not None:
self._dtype = self.data.dtype
return self._dtype
def __dtype_set(self, dtype):
if self.data is None:
self._dtype = dtype
else:
raise StateError('cannot set dtype after data has been set')
dtype = property(__dtype_get, __dtype_set)
# self._shape is used when self.data hasn't been set yet
def __shape_get(self):
if self.data is not None:
self._shape = self.data.shape
return self._shape
def __shape_set(self, shape):
if self.data is None:
self._shape = shape
else:
raise StateError('cannot set shape after data has been set')
shape = property(__shape_get, __shape_set)
def __add__(self, y): return add(self, y)
def __radd__(self, x): return add(x, self)
def __iadd__(self, y): return add_inplace(self, y)
def __sub__(self, y): return sub(self, y)
def __rsub__(self, x): return sub(x, self)
def __isub__(self, y): return sub_inplace(self, y)
def __mul__(self, y): return mul(self, y)
def __rmul__(self, x): return mul(x, self)
def __imul__(self, y): return mul_inplace(self, y)
def __div__(self, y): return div(self, y)
def __rdiv__(self, x): return div(x, self)
def __idiv__(self, y): return div_inplace(self, y)
def __pow__(self, y): return pow(self, y)
def __rpow__(self, x): return pow(x, self)
def __ipow__(self, y): return pow_inplace(self, y)
def __neg__(self): return neg(self)
T = property(lambda self: transpose(self))
Tc = property(lambda self: transpose_copy(self))
def __copy__(self): return array_copy(self)
def __getitem__(self, item): return get_slice(self, item)
def __getslice__(self, *args): return get_slice(self, slice(*args))
#################
# NumpyR Compatibility
#
spec = property(lambda self: (numpy.ndarray, self.dtype, self.shape))
def set_value_inplace(self, value):
if 0 == len(self.shape):
self.data.itemset(value) # for scalars
else:
self.data[:] = value # for matrices
self.state = gof.result.Computed
class _test_Numpy2(unittest.TestCase):
def setUp(self):
build_eval_mode()
numpy.random.seed(44)
def tearDown(self):
pop_mode()
def test_0(self):
r = Numpy2()
def test_1(self):
o = numpy.ones((3,3))
r = Numpy2(data=o)
self.failUnless(r.data is o)
self.failUnless(r.shape == (3,3))
self.failUnless(str(r.dtype) == 'float64')
def test_2(self):
r = Numpy2(data=[(3,3),'int32'])
self.failUnless(r.data is None)
self.failUnless(r.shape == (3,3))
self.failUnless(str(r.dtype) == 'int32')
r.alloc()
self.failUnless(isinstance(r.data, numpy.ndarray))
self.failUnless(r.shape == (3,3))
self.failUnless(str(r.dtype) == 'int32')
def test_3(self):
a = Numpy2(data=numpy.ones((2,2)))
b = Numpy2(data=numpy.ones((2,2)))
c = add(a,b)
self.failUnless(_approx_eq(c, numpy.ones((2,2))*2))
def test_4(self):
ones = numpy.ones((2,2))
a = Numpy2(data=ones)
o = numpy.asarray(a)
self.failUnless((ones == o).all())
def test_5(self):
ones = numpy.ones((2,2))
self.failUnless(_approx_eq(Numpy2(data=ones), Numpy2(data=ones)))
def cgen(name, behavior, names, vals, converters = None):
def cgetspecs(names, vals, converters):
d = {}
assert len(names) == len(vals)
for name, value in zip(names, vals):
d[name] = value.data
specs = weave.ext_tools.assign_variable_types(names, d, type_converters = converters) #, auto_downcast = 0)
return d, specs
if not converters:
converters = type_spec.default
for converter in converters:
assert isinstance(converter, type_spec.omega_type_converter_extension)
d, specs = cgetspecs(names, vals, converters)
template = {}
template['name'] = name
template['code'] = behavior
template['members'] = "".join([spec.struct_members_code() for spec in specs])
template['support'] = "".join([spec.struct_support_code() for spec in specs])
template['typedefs'] = "".join([spec.struct_typedefs() for spec in specs])
template['incref'] = "".join(["Py_INCREF(py_%s);\n" % spec.name for spec in specs if spec.use_ref_count])
template['decref'] = "".join(["Py_DECREF(py_%s);\n" % spec.name for spec in specs if spec.use_ref_count])
template['struct_contents'] = """
%(typedefs)s
%(members)s
%(support)s
void init(void) {
%(incref)s
}
void cleanup(void) {
%(decref)s
}
int execute(void) {
%(code)s
return 0;
}
""" % template
template['md5'] = md5.md5(template['struct_contents']).hexdigest()
template['struct_name'] = "_omega_%(name)s_%(md5)s" % template
struct = "struct %(struct_name)s { %(struct_contents)s\n};" % template
static = """
int %(struct_name)s_executor(%(struct_name)s* self) {
return self->execute();
}
void %(struct_name)s_destructor(void* executor, void* self) {
((%(struct_name)s*)self)->cleanup();
free(self);
}
""" % template
code = "%(struct_name)s* __STRUCT_P = new %(struct_name)s();\n" % template
code += "".join([spec.struct_import_code() for spec in specs])
code += "__STRUCT_P->init();\n"
code += "return_val = PyCObject_FromVoidPtrAndDesc((void*)(&%(struct_name)s_executor), __STRUCT_P, %(struct_name)s_destructor);\n" % template
return d, names, code, struct + static, converters
class Numpy2Op(gof.lib.PythonOp):
"""What can we do given we are interacting with Numpy2 inputs and outputs"""
def refresh(self, alloc = True):
shape = self.refresh_shape()
dtype = self.refresh_dtype()
out = self.out
if out.data is not None \
and out.shape == shape \
and out.dtype == dtype:
return
alloc |= out.data is not None
if alloc: out.data = None
out.shape = shape
out.dtype = dtype
if alloc: out.alloc()
class omega_op(Numpy2Op):
forbid_broadcast = False
@staticmethod
def __clsinit__(cls, name, bases, dct):
for fname in ['grad', 'c_impl', 'impl']:
if hasattr(cls, fname):
gof.make_static(cls, fname)
def __new__(cls, *inputs):
inputs = [wrap(input) for input in inputs]
return Numpy2Op.__new__(cls, *inputs)
def gen_outputs(self):
return [Numpy2() for i in xrange(self.nout)]
#TODO: use the version of this code that is in grad.py
# requires: eliminating module dependency cycles
def update_gradient(self, grad_d):
"""Call self.grad() and add the result to grad_d
This function is called by grad.Grad.bprop() to construct a symbolic gradient graph.
self.grad is called like this:
self.grad(*(self.inputs + [grad_d[output] for output in self.outputs]))
In general, grad() should return a list of ResultValue instances whose
length matches that of self.inputs, and whose elements are the
gradients of self.inputs.
There is a (but often used) special feature in place to automatically
wrap the return value of grad() in a list if it is a ResultValue instance
and the op is unary. This makes many grad implementations a little
cuter.
"""
inputgs = self.grad(*(self.inputs + [grad_d[output] for output in self.outputs]))
if len(self.inputs) == 1 and gof.result.is_result(inputgs):
inputgs = [inputgs]
else:
assert len(inputgs) == len(self.inputs)
for input, inputg in zip(self.inputs, inputgs):
grad_d.add(input, inputg)
def c_code(self, converters = None):
(inames, onames) = self.variable_names()
behavior = self._c_impl()
return cgen(self.__class__.__name__, behavior, inames + onames, self.inputs + self.outputs, converters)
def c_headers(self):
return []
def c_libs(self):
return []
def c_support_code(self):
return ""
def variable_names(self):
(inames, onames), _1, _2, _3 = inspect.getargspec(self.c_impl)
return (inames, onames)
def _c_impl(self):
return self.c_impl(self.inputs, self.outputs)
def c_impl(inputs, outputs):
raise NotImplementedError()
def c_compile_args(self):
# I always used these, but they don't make much improvement
#'-ffast-math', '-falign-loops=8'
return ['-O2']
def c_thunk_factory(self):
self.refresh()
d, names, code, struct, converters = self.c_code()
cthunk = object()
module_name = md5.md5(code).hexdigest()
mod = weave.ext_tools.ext_module(module_name)
instantiate = weave.ext_tools.ext_function('instantiate',
code,
names,
local_dict = d,
global_dict = {},
type_converters = converters)
instantiate.customize.add_support_code(self.c_support_code() + struct)
for arg in self.c_compile_args():
instantiate.customize.add_extra_compile_arg(arg)
for header in self.c_headers():
instantiate.customize.add_header(header)
for lib in self.c_libs():
instantiate.customize.add_library(lib)
#add_library_dir
#print dir(instantiate.customize)
#print instantiate.customize._library_dirs
if os.getenv('OMEGA_BLAS_LD_LIBRARY_PATH'):
instantiate.customize.add_library_dir(os.getenv('OMEGA_BLAS_LD_LIBRARY_PATH'))
mod.add_function(instantiate)
mod.compile(location = _compile_dir())
module = __import__("%s" % (module_name), {}, {}, [module_name])
def creator():
return module.instantiate(*[x.data for x in self.inputs + self.outputs])
return creator
def c_thunk(self):
return self.c_thunk_creator()
def c_perform(self):
thunk = self.c_thunk()
cutils.run_cthunk(thunk)
def upcast(dtype, *dtypes):
z = numpy.zeros((), dtype = dtype)
for dtype in dtypes:
z = z + numpy.zeros((), dtype = dtype)
return z.dtype
from grad import Undefined
def wrap_producer(f):
class producer(gof.lib.NewPythonOp):
def __init__(self, shape, dtype, order):
assert order == 'C' #TODO: let Numpy2 support this
if current_mode() == 'build_eval':
gof.lib.NewPythonOp.__init__(self,
[input(shape), input(dtype), input(order)],
[Numpy2(data = f(shape, dtype))])
elif current_mode() == 'build':
gof.lib.NewPythonOp.__init__(self,
[input(shape), input(dtype), input(order)],
[Numpy2(data = (shape, dtype))])
def gen_outputs(self):
return [Numpy2() for i in xrange(self.nout)]
impl = f
def grad(*args):
return [Undefined] * (len(args) - 1)
producer.__name__ = f.__name__
def ret(shape, dtype = 'float64', order = 'C'):
return producer(shape, dtype, order).out
return ret
ndarray = wrap_producer(numpy.ndarray)
array = wrap_producer(numpy.array)
zeros = wrap_producer(numpy.zeros)
ones = wrap_producer(numpy.ones)
class _testCase_producer_build_mode(unittest.TestCase):
def test_0(self):
"""producer in build mode"""
build_mode()
a = ones(4)
self.failUnless(a.data is None, a.data)
self.failUnless(a.state is gof.result.Empty, a.state)
self.failUnless(a.shape == 4, a.shape)
self.failUnless(str(a.dtype) == 'float64', a.dtype)
pop_mode()
def test_1(self):
"""producer in build_eval mode"""
build_eval_mode()
a = ones(4)
self.failUnless((a.data == numpy.ones(4)).all(), a.data)
self.failUnless(a.state is gof.result.Computed, a.state)
self.failUnless(a.shape == (4,), a.shape)
self.failUnless(str(a.dtype) == 'float64', a.dtype)
pop_mode()
if __name__ == '__main__':
unittest.main()
import unittest
from constructor import *
import random
class MyAllocator(Allocator):
def __init__(self, fn):
self.fn = fn
def __call__(self):
return self.fn.__name__
def f1(a, b, c):
return a + b + c
def f2(x):
return "!!%s" % x
class _test_Constructor(unittest.TestCase):
def test_0(self):
c = Constructor(MyAllocator)
c.update({"fifi": f1, "loulou": f2})
assert c.fifi() == 'f1' and c.loulou() == 'f2'
def test_1(self):
c = Constructor(MyAllocator)
c.add_module(random)
assert c.random.random() == 'random' and c.random.randint() == 'randint'
def test_2(self):
c = Constructor(MyAllocator)
c.update({"fifi": f1, "loulou": f2})
globals().update(c)
assert fifi() == 'f1' and loulou() == 'f2'
if __name__ == '__main__':
unittest.main()
from constructor import Allocator
from op import Op
class OpAllocator(Allocator):
def __init__(self, opclass):
if not issubclass(opclass, Op):
raise TypeError("Expected an Op instance.")
self.opclass = opclass
class FilteredOpAllocator(OpAllocator):
def filter(self, op):
pass
def __call__(self, *inputs):
op = self.opclass(*inputs)
self.filter(op)
if len(op.outputs) == 1:
return op.outputs[0]
else:
return op.outputs
class BuildAllocator(FilteredOpAllocator):
pass
class EvalAllocator(FilteredOpAllocator):
def filter(self, op):
op.perform()
for output in op.outputs:
output.role = None
class BuildEvalAllocator(FilteredOpAllocator):
def filter(self, op):
op.perform()
from utils import AbstractFunctionError
class Dispatcher(list):
all_dispatchers = {}
def __init__(self, name, description):
self.name = name
self.description = description
self.all_dispatchers[name] = self
def __call__(self, *inputs, **opts):
for candidate in self:
try:
return candidate(*inputs, **opts)
except TypeError:
continue
if opts:
s = " with options %s" % opts
else:
s = ""
raise OmegaTypeError("No candidate found for %s(%s) %s" \
% (self.name,
", ".join([input.__class__.__name__ for input in inputs]),
s))
def add_handler(self, x):
self.insert(0, x)
def fallback_handler(self, x):
self.append(x)
class Allocator:
def __init__(self, fn):
self.fn = fn
class IdentityAllocator(Allocator):
def __call__(self, *args, **kwargs):
return self.fn(*args, **kwargs)
class Constructor(dict):
def __init__(self, allocator):
self._allocator = allocator
def add_from_module(self, module):
for symbol in dir(module):
if symbol[:2] == '__': continue
obj = getattr(module, symbol)
try:
self[symbol] = self._allocator(obj)
except TypeError:
pass
def add_module(self, module, module_name = None):
if module_name is None:
module_name = module.__name__
d = Constructor(self._allocator)
d.add_from_module(module)
self[module_name] = d
def update(self, d, can_fail = False):
for name, fn in d.items():
self.add(name, fn, can_fail)
def add(self, name, fn, can_fail = True):
if isinstance(fn, Constructor):
self[name] = fn
else:
try:
self[name] = self._allocator(fn)
except TypeError:
if can_fail:
raise
def __getattr__(self, attr):
return self[attr]
from utils import OmegaError
class OmegaTypeError(OmegaError, TypeError):
pass
############################
# Dispatcher
############################
class Dispatcher(list):
all_dispatchers = {}
def __init__(self, name, description):
self.name = name
self.description = description
self.all_dispatchers[name] = self
def __call__(self, *inputs, **opts):
for candidate in self:
try:
return candidate(*inputs, **opts)
except OmegaTypeError:
continue
if opts:
s = " with options %s" % opts
else:
s = ""
raise OmegaTypeError("No candidate found for %s(%s) %s" \
% (self.name,
", ".join([input.__class__.__name__ for input in inputs]),
s))
def add_handler(self, x):
self.insert(0, x)
def fallback_handler(self, x):
self.append(x)
# Dispatchers for all python operators
Add = Dispatcher("Add", "x + y")
Subtract = Dispatcher("Subtract", "x - y")
Multiply = Dispatcher("Multiply", "x * y")
Divide = Dispatcher("Divide", "x / y")
FloorDivide = Dispatcher("FloorDivide", "x // y")
Modulo = Dispatcher("Modulo", "x % y")
Power = Dispatcher("Power", "x ** y")
Negate = Dispatcher("Negate", "-x")
Abs = Dispatcher("Abs", "abs(x)")
LeftShift = Dispatcher("LeftShift", "x << y")
RightShift = Dispatcher("RightShift", "x >> y")
Equals = Dispatcher("Equals", "x == y")
NotEquals = Dispatcher("NotEquals", "x != y")
Less = Dispatcher("Less", "x < y")
LessOrEqual = Dispatcher("LessOrEqual", "x <= y")
Greater = Dispatcher("Greater", "x > y")
GreaterOrEqual = Dispatcher("GreaterOrEqual", "x >= y")
Contains = Dispatcher("Contains", "x in y")
BinaryOr = Dispatcher("BinaryOr", "x | y")
BinaryAnd = Dispatcher("BinaryAnd", "x & y")
BinaryXor = Dispatcher("BinaryXor", "x ^ y")
BinaryInverse = Dispatcher("BinaryInverse", "~x")
# Dispatchers for special operations
Transpose = Dispatcher("Transpose", "x.T")
# Others
Log = Dispatcher("Log", 'log(x)')
Exp = Dispatcher("Exp", 'exp(x)')
Sin = Dispatcher("Sin", 'sin(x)')
Cos = Dispatcher("Cos", 'cos(x)')
Tan = Dispatcher("Tan", 'tan(x)')
############################
# PythonOperatorSupport
############################
class PythonOperatorSupport(object):
"""Support for built-in Python operators."""
# Common arithmetic operations
def __add__(self, x):
return Add(self, x)
def __radd__(self, x):
return Add(x, self)
def __sub__(self, x):
return Subtract(self, x)
def __rsub__(self, x):
return Subtract(x, self)
def __mul__(self, x):
return Multiply(self, x)
def __rmul__(self, x):
return Multiply(x, self)
def __div__(self, x):
return Divide(self, x)
def __rdiv__(self, x):
return Divide(x, self)
def __floordiv__(self, x):
return FloorDivide(self, x)
def __rfloordiv__(self, x):
return FloorDivide(x, self)
def __mod__(self, x):
return Modulo(self, x)
def __rmod__(self, x):
return Modulo(x, self)
def __pow__(self, x):
return Power(self, x)
def __rpow__(self, x):
return Power(x, self)
def __neg__(self):
return Negate(self)
def __abs__(self):
return Abs(self)
# Less common arithmetic operations
def __lshift__(self, x):
return LeftShift(self, x)
def __rlshift__(self, x):
return LeftShift(x, self)
def __rshift__(self, x):
return RightShift(self, x)
def __rrshift__(self, x):
return RightShift(x, self)
# Comparison operations
# def __eq__(self, x):
# return Equals(self, x)
# def __ne__(self, x):
# return NotEquals(self, x)
def __lt__(self, x):
return Less(self, x)
def __le__(self, x):
return LessOrEqual(self, x)
def __gt__(self, x):
return Greater(self, x)
def __ge__(self, x):
return GreaterOrEqual(self, x)
def __contains__(self, x):
return Contains(self, x)
# Binary operations
def __or__(self, x):
return BinaryOr(self, x)
def __ror__(self, x):
return BinaryOr(x, self)
def __and__(self, x):
return BinaryAnd(self, x)
def __rand__(self, x):
return BinaryAnd(x, self)
def __xor__(self, x):
return BinaryXor(self, x)
def __rxor__(self, x):
return BinaryXor(x, self)
def __invert__(self):
return BinaryInverse(self)
# Other operations
T = property(lambda self: Transpose(self))
norm = property(lambda self: Norm(self))
# Always nonzero
def __nonzero__(self):
return True
__all__ = globals().keys()
......@@ -110,9 +110,9 @@ class Env(graph.Graph):
### Public interface ###
# def add_output(self, output):
# self.outputs.add(output)
# self.__import_r__([output])
def add_output(self, output):
self.outputs.add(output)
self.__import_r__([output])
def clients(self, r):
"Set of all the (op, i) pairs such that op.inputs[i] is r."
......@@ -250,8 +250,6 @@ class Env(graph.Graph):
if r in self.outputs:
was_output = True
self.outputs[self.outputs.index(r)] = new_r
# self.outputs.remove(r)
# self.outputs.add(new_r)
# The actual replacement operation occurs here. This might raise
# an error.
......@@ -263,8 +261,6 @@ class Env(graph.Graph):
if was_output:
if not new_was_output:
self.outputs[self.outputs.index(new_r)] = r
# self.outputs.remove(new_r)
# self.outputs.add(r)
# Move back the clients. This should never raise an error.
self.__move_clients__(clients, new_r, r)
......
......@@ -118,7 +118,6 @@ class DestroyHandler(Listener, Constraint, Orderings):
for input in destroyed:
path = self.__path__(input)
self.__add_destroyer__(path + [output])
####### self.__add_destroyer__(path + [op])
elif views:
if len(views) > 1:
......
from copy import copy
from op import Op
from result import is_result, ResultBase
from utils import ClsInit, Keyword, AbstractFunctionError
import opt
import env
import features
import ext
from python25 import all
__all__ = [ 'UNDEFINED',
'current_mode',
'set_mode',
'build_mode',
'eval_mode',
'build_eval_mode',
'pop_mode',
'DummyOp',
'DummyRemover',
'PythonOp',
'PythonOpt',
'make_static']
UNDEFINED = Keyword("UNDEFINED", False)
def make_static(cls, fname):
f = getattr(cls, fname)
if hasattr(f, 'im_func'):
f = f.im_func
setattr(cls, fname, staticmethod(f))
def compute_from(nodes, history):
"""Recursively evaluate each node (in a quick & dirty way).
history (aka inputs) is a set of nodes that need not be [re]computed.
TODO: make this more correct by building a little graph and executing it.
The current implementation doesn't take into account any ordering
constraints imposed by destructors, for example.
"""
def compute_recursive(node):
if node and (node not in history):
if hasattr(node, 'owner'): #node is storage
compute_recursive(node.owner)
else: #node is op
if node.destroy_map():
raise ValueError('compute_from() does not work on nodes with destroy_maps')
for input in node.inputs:
compute_recursive(input)
node.perform()
history.add(node)
for n in nodes:
compute_recursive(n)
def compute(*nodes):
"""Recursively evaluate each node (in a quick & dirty way)."""
compute_from(nodes, set())
def root_inputs(input):
"""Return the leaves of a search through consecutive view_map()s"""
owner = input.owner
if owner:
view_map = owner.view_map()
if input in view_map:
answer = []
for input2 in view_map[input]:
answer.extend(root_inputs(input2))
return answer
else:
return [input]
else:
return [input]
class ForbidConstantOverwrite(features.Listener, features.Constraint):
def __init__(self, env):
self.env = env
self.bad = set()
def on_import(self, op):
for output, inputs in op.destroy_map().items():
for input in inputs:
for root_input in root_inputs(input):
if getattr(root_input, 'constant', False):
self.bad.add(op)
return
def on_prune(self, op):
if op in self.bad:
self.bad.remove(op)
def on_rewire(self, clients, r, new_r):
for op, i in clients:
self.on_prune(op)
self.on_import(op)
def validate(self):
if self.bad:
raise env.InconsistencyError("The following ops overwrite a constant value: %s" % self.bad)
else:
return True
class NewPythonOp(Op):
__env_require__ = DestroyHandler, ForbidConstantOverwrite
def view_map(self):
return {}
def destroy_map(self):
return {}
class PythonOp(NewPythonOp):
__metaclass__ = ClsInit
nout = 1
@staticmethod
def __clsinit__(cls, name, bases, dct):
# make impl a static method
cls.set_impl(cls.impl)
def __new__(cls, *inputs, **kwargs):
op = NewPythonOp.__new__(cls)
op.__init__(*inputs)
mode = kwargs.get('mode', None) or current_mode()
if mode == 'eval':
op.perform()
if op.nout == 1:
return op.out.data
else:
return [output.data for output in op.outputs]
elif mode == 'build_eval':
op.perform()
if op.nout == 1:
return op.out
else:
return op.outputs
def __init__(self, *inputs):
NewPythonOp.__init__(self, inputs, self.gen_outputs())
def __validate__(self):
return all([is_result(i) for i in self.inputs])
def gen_outputs(self):
raise AbstractFunctionError()
def check_input(self, input):
def input_is_up_to_date(input):
answer = True
for input in root_inputs(input):
answer &= input.up_to_date
return answer
if input.data is None:
raise ValueError("Uncomputed input: %s in %s" % (input, self))
if not input_is_up_to_date(input):
raise ValueError("Input is out of date: %s in %s" % (input, self))
def perform(self):
def input_is_constant(input):
answer = False
for input in root_inputs(input):
answer |= input.constant
return answer
exc = set()
for output, inputs in self.destroy_map().items():
exc.update(inputs)
for input in inputs:
if self.input_is_constant(input):
raise ValueError("Input is constant: %s" % input)
for input in exc:
self.check_input(input)
input.up_to_date = False
for input in self.inputs:
if input not in exc:
self.check_input(input)
if 0:
#J- why is this try catch here? Leftover debug?
try:
results = self._impl()
except Exception, e:
print "Error in %s: %s" % (self, e)
raise
else:
results = self._impl()
if self.nout == 1:
self.out.data = results
else:
assert self.nout == len(results)
for result, output in zip(results, self.outputs):
output.data = result
def _perform(self):
results = self._impl()
if self.nout == 1:
self.out.data = results
else:
assert self.nout == len(results)
for result, output in zip(results, self.outputs):
output.data = result
def _perform_inplace(self):
results = self._impl()
if self.nout == 1:
self.out.set_value_inplace(results)
else:
assert self.nout == len(results)
for result, output in zip(results, self.outputs):
output.set_value_inplace(result)
def _impl(self):
return self.impl(*[input.data for input in self.inputs])
@classmethod
def set_impl(cls, impl):
make_static(cls, 'impl')
def impl(*args):
raise NotImplementedError("This op has no implementation.")
def __copy__(self):
"""
Copies the inputs list shallowly and copies all the outputs
because of the one owner per output restriction.
"""
# new_inputs = copy(op.inputs)
# # We copy the outputs because they are tied to a single Op.
# new_outputs = [copy(output) for output in op.outputs]
build_mode()
op = self.__class__(*self.inputs)
pop_mode()
# op._inputs = new_inputs
# op._outputs = new_outputs
# for i, output in enumerate(op.outputs):
# # We adjust _owner and _index manually since the copies
# # point to the previous op (self).
# output._owner = op
# output._index = i
if isinstance(op, (list, tuple)):
return op[0].owner
return op.owner
__mode__ = ['build_eval']
def current_mode():
return __mode__[-1]
def set_mode(mode):
__mode__.append(mode)
def build_mode():
set_mode('build')
def eval_mode():
set_mode('eval')
def build_eval_mode():
set_mode('build_eval')
def pop_mode():
if len(__mode__) == 1:
raise Exception("There's only one mode left on the stack.")
else:
__mode__.pop()
class PythonOpt(opt.Optimizer):
def __init__(self, opt):
self.opt = opt
def optimize(self, env):
build_mode()
self.opt.optimize(env)
pop_mode()
class DummyOp(NewPythonOp):
def __init__(self, input):
Op.__init__(self, [input], [ResultBase()])
def thunk(self):
return lambda:None
DummyRemover = opt.OpRemover(DummyOp)
if 0:
class RefreshableOp(NewPythonOp):
def _specs(self):
try:
return self.specs(*[input.spec for input in self.inputs])
except NotImplementedError:
raise NotImplementedError("%s cannot infer the specs of its outputs" % self.__class__.__name__)
def specs(*inputs):
raise NotImplementedError
def refresh(self):
"""Update and allocate outputs if necessary"""
for input in self.inputs:
input.refresh()
change = self._propagate_specs()
if change:
self.alloc(except_list)
return change
def _propagate_specs(self):
specs = self._specs()
if self.nout == 1:
specs = [specs]
change = False
for output, spec in zip(self.outputs, specs):
if output.spec != spec:
output.spec = spec
change = True
return change
def alloc(self, except_list = []):
for output in self.outputs:
if output not in except_list:
output.alloc()
......@@ -7,10 +7,6 @@ __all__ = ['ModalConstructor',
'build',
'eval',
'build_eval',
# 'ModalWrapper',
# 'BuildMode',
# 'EvalMode',
# 'BuildEvalMode',
'make_constructors',
]
......@@ -31,28 +27,15 @@ class ModalConstructor:
if mode != modal_wrapper:
raise TypeError("Inconsistent modes.")
fn_args.append(arg)
# for arg in args:
# if isinstance(arg, ModalWrapper):
# if modal_wrapper is None:
# modal_wrapper = arg.__class__
# else:
# if not isinstance(arg, modal_wrapper):
# raise TypeError("Inconsistent modes.")
# fn_args.append(arg.r)
# else:
# fn_args.append(arg)
op = self.fn(*fn_args)
if modal_wrapper:
modal_wrapper(op)
# modal_wrapper.filter(op)
for output in op.outputs:
output.__mode__ = modal_wrapper
if len(op.outputs) == 1:
return op.outputs[0]
#return modal_wrapper(op.outputs[0])
else:
return op.outputs
#return [modal_wrapper(output) for output in op.outputs]
def add_modal_members(cls, *members):
......@@ -65,30 +48,6 @@ def add_modal_members(cls, *members):
setattr(cls, member, fn(member))
# class ModalWrapper:
# def __init__(self, r):
# self.r = r
# def __as_result__(self):
# return self.r
# def __get_owner(self):
# return self.r.owner
# owner = property(__get_owner)
# @classmethod
# def filter(cls, op):
# raise AbstractFunctionError()
# members1 = 'add sub mul div pow floordiv mod pow lshift rshift and or xor'.split(' ')
# members = []
# members += ["__%s__" % x for x in members1 + 'neg invert'.split(' ')]
# members += ["__r%s__" % x for x in members1]
# add_modal_members(ModalWrapper, *members)
def build_mode(op):
pass
......@@ -112,24 +71,6 @@ eval = mode_setter(eval_mode)
build_eval = mode_setter(build_eval_mode)
# class BuildMode(ModalWrapper):
# @classmethod
# def filter(cls, op):
# pass
# class EvalMode(ModalWrapper):
# @classmethod
# def filter(cls, op):
# op.perform()
# for output in op.outputs:
# output._role = None
# class BuildEvalMode(ModalWrapper):
# @classmethod
# def filter(cls, op):
# op.perform()
def _is_op(x):
try: return issubclass(x, Op)
except: return False
......
......@@ -142,22 +142,6 @@ class Op(object):
"""
raise AbstractFunctionError()
# def c_update(self):
# """
# Returns C code that allocates and/or updates the outputs
# (eg resizing, etc.) so they can be manipulated safely
# by c_code.
# You may use the variable names defined by c_var_names()
# """
# raise AbstractFunctionError()
# def c_update_cleanup(self):
# """
# Clean up things allocated by c_update().
# """
# raise AbstractFunctionError()
def c_code(self):
"""
Returns C code that does the computation associated to this
......
from env import Env
from utils import AbstractFunctionError
class Prog:
def __init__(self, inputs, outputs, optimizer, linker_class, features = []):
self.inputs = inputs
if isinstance(outputs, dict):
for name, output in outputs.items():
setattr(self, name, output)
self.outputs = outputs.values()
else:
self.outputs = outputs
self.optimizer = optimizer
self.env = Env(self.inputs, self.outputs, features, False)
self.env.add_feature(EquivTool)
self.linker = linker_class(self.env)
def build(self):
self.optimizer.optimize(self.env)
def equiv(self, r):
return self.env.equiv(r)
def __getitem__(self, r):
if isinstance(r, str):
return getattr(self, r)
else:
return self.equiv(r)
def __setitem__(self, r, value):
if isinstance(r, tuple):
for a, b in zip(r, value):
self.__setitem__(a, b)
else:
self[r].data = value
# import compile
import env
import link
from features import EquivTool
class Prog:
def __init__(self, inputs, outputs, optimizer, linker, features = []):
self.optimizer = optimizer
self.linker = linker
features = set(features)
features.add(EquivTool)
self.env = env.Env(inputs, outputs, features) #, False)
self.optimizer.optimize(self.env)
self.perform = self.linker(self.env)
self.outputs = outputs
# def __optimize__(self):
# self.optimizer.apply(self.env)
# self.order = self.env.toposort()
def equiv(self, r):
return self.env.equiv(r)
def __getitem__(self, r):
return self.equiv(r)
def __setitem__(self, r, value):
if isinstance(r, tuple):
for a, b in zip(r, value):
self.__setitem__(a, b)
else:
self.equiv(r).set_value(value)
def __call__(self, *args):
self.perform()
for output in self.outputs:
output.set_value(self[output])
return self.outputs
# return [output for output in self.env.outputs]
# if args:
# for input, arg in zip(self.inputs, args):
# if arg is not None:
# input.value = arg
# for thunk, op in zip(self.thunks, self.order):
# try:
# thunk()
# except Exception, e:
# raise e.__class__("Error in " + str(op) + ": " + str(e))
# return [output.value for output in self.outputs]
import sys
if sys.version_info[:2] < (2,5):
def all(iterable):
for element in iterable:
if not element:
return False
return True
else:
# Only bother with this else clause and the __all__ line if you are putting
# this in a separate file.
import __builtin__
all = __builtin__.all
__all__ = ['all']
......@@ -160,32 +160,17 @@ class ResultBase(object):
"""
raise AbstractFunctionError()
# def c_extract(self):
# get_from_list = """
# //PyObject* py_%(name)s = PyList_GET_ITEM(%(name)s_storage, 0);
# //Py_XINCREF(py_%(name)s);
# """
# return get_from_list + self.c_data_extract()
def c_extract(self):
"""
# The code returned from this function must be templated using
# "%(name)s", representing the name that the caller wants to
# call this Result. The Python object self.data is in a
# variable called "py_%(name)s" and this code must declare a
# variable named "%(name)s" of type "%(type)s" where "%(type)s"
# will be replaced by the return value of
# self.c_type(). Additional variables and typedefs may not be
# produced. If the data is improper, set an appropriate error
# message and insert "%(fail)s".
The code returned from this function must be templated using
"%(name)s", representing the name that the caller wants to
call this Result. The Python object self.data is in a
variable called "py_%(name)s" and this code must set the
variables declared by c_declare to something representative
of py_%(name)s. If the data is improper, set an appropriate error
message and insert "%(fail)s".
"""
raise AbstractFunctionError()
# def c_cleanup(self):
# decref = """
# //Py_XDECREF(py_%(name)s);
# """
# return self.c_data_cleanup() + decref
def c_cleanup(self):
"""
......@@ -196,13 +181,6 @@ class ResultBase(object):
Note: EITHER c_cleanup OR c_sync will be called.
"""
raise AbstractFunctionError()
# def c_sync(self):
# set_in_list = """
# //PyList_SET_ITEM(%(name)s_storage, 0, py_%(name)s);
# //Py_XDECREF(py_%(name)s);
# """
# return self.c_data_sync() + set_in_list
def c_sync(self):
"""
......@@ -287,38 +265,3 @@ class ResultBase(object):
def same_properties(self, other):
raise AbstractFunction()
#################
# NumpyR Compatibility
#
up_to_date = property(lambda self: True)
def refresh(self): pass
def set_owner(self, owner, idx):
self.role = (owner, idx)
def set_value(self, value):
self.data = value #may raise exception
# def c_data_extract(self):
# return """
# PyArrayObject* %%(name)s;
# if (py_%%(name)s == Py_None)
# %%(name)s = NULL;
# else
# %%(name)s = (PyArrayObject*)(py_%%(name)s);
# typedef %(dtype)s %%(name)s_dtype;
# """ % dict(dtype = self.dtype)
# def c_data_sync(self):
# return """
# if (!%(name)s) {
# Py_XDECREF(py_%(name));
# py_%(name)s = Py_None;
# }
# else if ((void*)py_%(name)s != (void*)%(name)s) {
# Py_XDECREF(py_%(name));
# py_%(name)s = (PyObject*)%(name)s;
# }
# """
#ifndef _OMEGA_H
#define _OMEGA_H
//#include whatever defines PyArrayObject
template<typename T>
struct TMat_t
{
T * __restrict__ d;/**< pointer to element (0,0) */
size_t M; /**< number of rows */
size_t N; /**< number of columns */
size_t m; /**< row stride */
size_t n; /**< column stride */
bool invalid;
/** null */
TMat_t(const PyArrayObject *o) :
d((double*) o->data),
M((o->nd==2) ? o->dimensions[0] : 0),
N((o->nd==2) ? o->dimensions[1] : 0),
m((o->nd==2) ? o->strides[0] / sizeof(double) : 0),
n((o->nd==2) ? o->strides[1] / sizeof(double) : 0),
invalid((o->nd !=2) || (o->descr->elsize != sizeof(T)))
{
}
/** unsafe element access */
const T & operator()(size_t i, size_t j) const
{
return d[ i * m + j*n];
}
/** unsafe element access */
T & operator()(size_t i, size_t j)
{
return d[ i * m + j*n];
}
/** safe element access */
const T & at(size_t i, size_t j) const
{
return d[ assert((i < M) && (j < N)), i * m + j*n];
}
/** safe element access */
T & at(size_t i, size_t j)
{
return d[ assert((i < M) && (j < N)), i * m + j*n];
}
};
#endif
from elemwise import elemwise
def scalar_switch(normal_f, scalar_f, scalar_f_reverse = None):
def f(x, y):
x, y = wrap(x), wrap(y)
if y.constant and not y.data.shape:
return scalar_f(x, y)
if x.constant and not x.data.shape:
if scalar_f_reverse:
return scalar_f_reverse(y, x)
else:
raise TypeError("You cannot do this operation on a scalar.")
return normal_f(x, y)
return f
# Wrapper to ensure that all inputs to the function impl have the same size (foils numpy's broadcasting)
def assert_same_shapes(impl):
def ret(x, *rest):
shape = x.shape
for other in rest:
if other.shape != shape:
raise ValueError("The dimensions of the inputs do not match.")
return impl(x, *rest)
return ret
# Wrapper to ensure that the last input to impl is a scalar
def tensor_scalar_impl(impl):
def ret(x, a):
if a.shape:
raise ValueError("The second argument to %s must be a scalar." % impl)
return impl(x, a)
return ret
class tensor_scalar_op(elemwise):
@classmethod
def variable_names(cls):
return (['x', '_a'], ['z', ])
@classmethod
def loop_variables(cls):
return (['x', ], ['z', ])
def c_init((x, _a), (z, )):
return """
if (PyArray_SIZE(_a) != 1) {
PyErr_SetString(PyExc_ValueError, \"The size of the scalar argument is not 1.\");
}
_a_dtype a = ((_a_dtype*)PyArray_DATA(_a))[0];
"""
def _c_foreach(self):
return "z_i = %s;" % self.c_expr
## Addition ##
class add_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__add__)
def grad(x, y, gz):
return gz, gz
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i + y_i;"
add_elemwise_inplace = add_elemwise.inplace_version()
add_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__iadd__))
class add_scalar(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__add__)
def grad(x, a, gz):
return gz, sum(gz)
c_expr = "x_i + a"
add_scalar_inplace = add_scalar.inplace_version()
add_scalar_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__iadd__))
class twice(elemwise):
def impl(x):
return 2.0 * x
def grad(x, gz):
return scale(gz, 2.0)
def c_foreach((x_i, ), (z_i, )):
"z_i = x_i + x_i;"
twice_inplace = twice.inplace_version()
## Subtraction ##
class sub_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__sub__)
def grad(x, y, gz):
return gz, -gz
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i - y_i;"
sub_elemwise_inplace = sub_elemwise.inplace_version()
sub_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__isub__))
def sub_scalar_r(x, a):
return add_scalar(x, -a)
def sub_scalar_l(x, a):
return add_scalar(-x, a)
def sub_scalar_r_inplace(x, a):
return add_scalar_inplace(x, -a)
## Element-wise multiplication ##
class mul_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__mul__)
def grad(x, y, gz):
return mul(y, gz), mul(x, gz)
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i * y_i;"
mul_elemwise_inplace = mul_elemwise.inplace_version()
mul_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__imul__))
class scale(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__mul__)
def grad(x, a, gz):
return scale(a, gz), sum(mul_elemwise(x, gz))
c_expr = "x_i * a"
scale_inplace = scale.inplace_version()
scale_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__imul__))
class sqr(elemwise):
def impl(x):
return x * x
def grad(x, gz):
return scale(mul_elemwise(x, gz), 2.0)
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i * x_i;"
isqr = sqr.inplace_version()
isqr.set_impl(lambda x: x.__imul__(x))
class sqrt(elemwise):
impl = numpy.sqrt
def grad(x, gz):
return scale(div(gz, sqrt(x)), 0.5)
def c_foreach((x_i, ), (z_i, )):
return "z_i = pow(x_i, 0.5);"
isqrt = sqrt.inplace_version()
isqrt.set_impl(lambda x: x.__ipow__(0.5))
## Element-wise division ##
class div_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__div__)
def grad(x, y, gz):
return div(gz, y), -div(mul(x, gz), sqr(y))
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i / y_i;"
div_elemwise_inplace = div_elemwise.inplace_version()
div_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__idiv__))
def div_scalar_r(x, a):
return scale(x, inv_elemwise(a))
def div_scalar_l(x, a):
return scale(inv_elemwise(x), a)
def div_scalar_r_inplace(x, a):
return scale_inplace(x, inv_elemwise(a))
## Scaling ##
class scale(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__mul__)
def grad(x, a, gz):
return scale(a, gz), sum(mul_elemwise(x, gz))
c_expr = "x_i * a"
scale_inplace = scale.inplace_version()
scale_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__imul__))
class neg(elemwise):
impl = numpy.ndarray.__neg__
def grad(x, gz):
return -gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = -x_i;"
neg_inplace = neg.inplace_version()
neg_inplace.set_impl(lambda x: x.__imul__(-1))
class inv_elemwise(elemwise):
impl = lambda x: 1 / x
def grad(x, gz):
return -gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = 1 / x_i;"
inv_elemwise_inplace = inv_elemwise.inplace_version()
## Dot product ##
class dot(omega_op):
@staticmethod
def _output_shape(xshape, yshape):
if len(xshape) == 0: # x is a scalar
shape = yshape
else:
if len(yshape) >= 2: #y is a matrix or tensor
assert xshape[-1] == yshape[-2]
shape = tuple(xshape[:-1]+ yshape[:-2]+yshape[-1:])
elif len(yshape)==1: #y is vector
assert xshape[-1] == yshape[-1]
shape = tuple(xshape[:-1])
else: #y is a scalar
shape = xshape
return shape
impl = numpy.dot
def grad(x, y, gz):
return dot(gz, transpose(y)), dot(transpose(x), gz)
def refresh(self, alloc=False):
x,y = self.inputs
shape = self._output_shape(x.shape, y.shape)
dtype = upcast(x.dtype, y.dtype)
if self.out.data is not None \
and self.out.shape == shape \
and self.out.dtype == dtype:
return #everything is ok
if alloc or self.out.data is not None: #data should be allocated
self.out.data = None
self.out.shape = shape
self.out.dtype = dtype
self.out.alloc()
else:
self.out.shape = shape
self.out.dtype = dtype
def c_support_code(self):
return blas.cblas_header_text()
def c_libs(self):
return blas.ldflags()
def c_impl((_x, _y), (_z, )):
return blas.gemm_code('', '1.0', '0.0')
## Transposition ##
class transpose(omega_op):
def view_map(self): return {self.out: [self.inputs[0]]}
impl = numpy.transpose
def grad(x, gz):
return transpose_copy(gz)
def refresh_shape(self):
rval = list(self.inputs[0].shape)
rval.reverse()
return rval
def refresh_dtype(self):
return self.inputs[0].dtype
def c_impl((x, ), (xt, )):
return """
const int l = x->nd;
// The user must ensure that all references to
//xt->data go through xt, or there's going to be trouble..
int refcheck = 0;
if (x == xt)
{
return -1;
}
if (refcheck)
{
int refcnt = PyArray_REFCOUNT(xt);
if ((refcnt > 2) // you might think this should be 1.. but this works
//|| (xt->base != NULL)
|| (xt->weakreflist != NULL))
{
PyErr_SetString(PyExc_ValueError,
"cannot resize an array that has "\\
"been referenced or is referencing\\n"\\
"another array in this way. Use the "\\
"resize function");
return -2;
}
}
if (xt->nd != x->nd)
{
// this technique comes from PyArray_Resize()
npy_intp * dimptr = (npy_intp*)PyDimMem_RENEW(xt->dimensions, 2 * x->nd);
if (!dimptr)
{
PyErr_NoMemory();
return 1;
}
xt->nd = x->nd;
xt->dimensions = dimptr;
xt->strides = dimptr + x->nd;
}
//copy x's dimensions and strides
for (int i = 0; i < l; ++i)
{
xt->dimensions[i] = x->dimensions[l-i-1];
xt->strides[i] = x->strides[l-i-1];
}
// point directly at b's type descriptor
Py_INCREF(x->descr);
Py_DECREF(xt->descr);
xt->descr = x->descr;
// name x as a base of xt, increment its refcount
if ( xt->base != (PyObject*)x)
{
Py_INCREF(x);
if ((xt->base) && (xt->base != Py_None))
{
Py_DECREF(xt->base);
}
xt->base = (PyObject*)x;
}
// mark xt as not owning its data
if (PyArray_CHKFLAGS(xt, NPY_OWNDATA))
{
PyDataMem_FREE(xt->data);
xt->flags &= ~NPY_OWNDATA;
}
xt->data = x->data;
// this function is described in
// ~/zzz.NOBACKUP/pub/src/numpy-1.0.3.1/numpy/core/src/arrayobject.c:1890
PyArray_UpdateFlags(xt, NPY_CONTIGUOUS|NPY_FORTRAN|NPY_ALIGNED|NPY_WRITEABLE);
/*
TODO
What should be done with the weakreflist ?
*/
"""
def transpose_copy(x):
return array_copy(transpose(x))
## Copy ##
class array_copy(elemwise):
impl = numpy.array
grad = lambda x, gz: gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i;"
## Power ##
class sqr(elemwise):
def impl(x):
return x * x
def grad(x, gz):
return scale(mul_elemwise(x, gz), 2.0)
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i * x_i;"
sqr_inplace = sqr.inplace_version()
sqr_inplace.set_impl(lambda x: x.__imul__(x))
class sqrt(elemwise):
impl = numpy.sqrt
def grad(x, gz):
return scale(div(gz, sqrt(x)), 0.5)
def c_foreach((x_i, ), (z_i, )):
return "z_i = pow(x_i, 0.5);"
sqrt_inplace = sqrt.inplace_version()
sqrt_inplace.set_impl(lambda x: x.__ipow__(0.5))
class exp(elemwise):
def impl(x): return numpy.exp(x)
def grad(x, gz): return gz * exp(x)
def c_foreach((x_i, ), (z_i, )): return "z_i = exp(x_i);"
class log(elemwise):
def impl(x): return numpy.log(x)
def grad(x, gz): return gz / x
def c_foreach((x_i, ), (z_i, )): return "z_i = log(x_i);"
class log2(elemwise):
def impl(x): return numpy.log2(x)
def grad(x, gz): return gz / (x * numpy.log(2))
def c_foreach((x_i, ), (z_i, )): return "z_i = log2(x_i);"
class pow_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__pow__)
def grad(x, s, gz):
raise NotImplemented # no gs
return gz * s * (pow_elemwise(x, s-1.0))
def c_foreach((x_i, s_i), (z_i, )):
return "z_i = pow(x_i, s_i)"
pow_elemwise_inplace = pow_elemwise.inplace_version()
pow_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__ipow__))
class pow_scalar_l(tensor_scalar_op):
impl = tensor_scalar_impl(lambda x, y: numpy.ndarray.__pow__(y, x))
def grad(x, s, gz):
raise NotImplemented # no gs
return gz * x * (pow_scalar_l(s,x-1.0))
c_expr = "pow(a, x_i)"
class pow_scalar_r(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__pow__)
def grad(x, s, gz):
gx = gz * s * (pow_scalar_r(x,s-1.0))
gs = sum(gz * pow_scalar_r(x,s) * log(x))
return gx, gs
c_expr = "pow(x_i, a)"
pow_scalar_r_inplace = pow_scalar_r.inplace_version()
pow_scalar_r_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__ipow__))
## Others ##
class minmax(elemwise):
nout = 2
def impl(x):
return x.min, x.max
def specs(x):
return [(numpy.ndarray, x[1], ())] * 2
# def alloc((x, ), (_min, _max)):
# _min.data = numpy.ndarray((), x.dtype)
# _max.data = numpy.ndarray((), x.dtype)
def c_init((x, ), (_min, _max)):
raise NotImplementedError
return """
_x_dtype min = _x[0];
_x_dtype max = _x[0];
"""
def c_foreach((x, ), (_min, _max)):
return """
if (x < min) min = x;
if (x > max) max = x;
"""
def c_finalize((x, ), (_min, _max)):
return """
_min[0] = min;
_max[0] = max;
"""
class fill(elemwise):
impl = lambda model, value: (model * 0) + value
def c_init((model, value), (z, )):
return "value_dtype value0 = ((value_dtype*)PyArray_DATA(value))[0];"
def c_foreach((model_i, value), (z_i, )):
return "z_i = value0;"
fill_inplace = fill.inplace_version()
class sum(elemwise):
impl = numpy.sum
def grad(x, gz):
return fill(x, gz)
def refresh_shape(self):
return ()
def c_init((x, ), (sum, )):
return "sum_dtype* sump = ((sum_dtype*)PyArray_DATA(sum)); sump[0] = 0;"
def c_foreach((x_i, ), (sum, )):
return "sump[0] += x_i;"
class ones_like(elemwise):
impl = numpy.ones_like
def grad(x, gz): return Undefined
class zeros_like(elemwise):
impl = numpy.zeros_like
def grad(x, gz): return Undefined
## Array slicing ##
class get_slice(omega_op):
def view_map(self): return {self.out: [self.inputs[0]]}
def impl(x, item):
rval = x.__getitem__(item)
#print 'get_slice running', rval
return rval
def grad(x, gz): raise NotImplemented
def refresh_shape(self):
x,item = self.inputs
rval = x.data.__getitem__(item.data).shape
#print 'refresh_shape', rval
return rval
def refresh_dtype(self):
return self.inputs[0].data.dtype
add = scalar_switch(add_elemwise, add_scalar, add_scalar)
add_inplace = scalar_switch(add_elemwise_inplace, add_scalar_inplace)
sub = scalar_switch(sub_elemwise, sub_scalar_r, sub_scalar_l)
sub_inplace = scalar_switch(sub_elemwise_inplace, sub_scalar_r_inplace)
mul = scalar_switch(mul_elemwise, scale, scale)
mul_inplace = scalar_switch(mul_elemwise_inplace, scale_inplace)
div = scalar_switch(div_elemwise, div_scalar_r, div_scalar_l)
div_inplace = scalar_switch(div_elemwise_inplace, div_scalar_r_inplace)
pow = scalar_switch(pow_elemwise, pow_scalar_r, pow_scalar_l)
pow_inplace = scalar_switch(pow_elemwise_inplace, pow_scalar_r_inplace)
import core
import gof
from numpy import random as r
# def rwrap(f):
# wrapped =
# def ret(self, *args):
class RandomState(gof.Op, gof.ext.IONames):
input_names = ['seed']
def __init__(self, seed):
inputs = [wrap(seed)]
outputs = [ResultValue()]
gof.Op.__init__(self, inputs, outputs)
def thunk(self):
def f():
self.out.storage = r.RandomState(self.seed.storage)
return f
class Random(object):
def __init__(seed):
self.state = core.wrap(seed)
......@@ -23,11 +23,30 @@ col = _broadcastable_pattern([0, 1])
class Tensor(ResultBase):
def __init__(self, dtype, broadcastable, data=None, name=None):
def __init__(self, dtype=None, broadcastable=None, data=None, name=None, constant=False):
if dtype is None or broadcastable is None:
if data is None:
raise TypeError("Provide non-None data to complete the dtype and broadcastable flags.")
dtype = data.dtype
if constant:
broadcastable = [1*(x == 1) for x in data.shape]
else:
broadcastable = [0] * len(data.shape)
self.broadcastable = broadcastable
self.dtype = str(dtype)
self.constant = constant
ResultBase.__init__(self, role = None, data = None, name = name)
def __get_constant(self):
return self._constant
def __set_constant(self, value):
if value:
self.indestructible = True
self._constant = value
constant = property(__get_constant, __set_constant)
def filter(self, data):
arr = numpy.asarray(data, dtype = self.dtype)
for b, s in zip(self.broadcastable, arr.shape):
......
......@@ -9,6 +9,12 @@ def upcast(dtype, *dtypes):
z = z + numpy.zeros((), dtype = dtype)
return str(z.dtype)
def wrap_as_tensor(x):
if isinstance(x, Tensor):
return x
else:
return Tensor(data=x, constant=True)
class TensorOp(Op):
nin = -1
......@@ -18,12 +24,6 @@ class TensorOp(Op):
def __init__(self, *inputs):
def wrap_as_tensor(x):
if isinstance(x, Tensor):
return x
else:
return Tensor(x)
inputs = map(wrap_as_tensor, inputs)
if self.nin >= 0:
......@@ -76,7 +76,6 @@ class BinaryTensorOp(TensorOp):
nin = 2
class Transpose(UnaryTensorOp):
def propagate_broadcastable(self, x):
......@@ -90,16 +89,562 @@ class Transpose(UnaryTensorOp):
def c_impl(self, x, z):
return """
PyArrayObject* transposed = (PyArrayObject*)PyArray_Transpose(%(x)s, NULL);
if (PyArray_REFCOUNT(transposed) == 1) {
printf("lala\\n");
}
if (%(z)s) {
Py_XDECREF(%(z)s);
}
//if (PyArray_REFCOUNT(transposed) == 1) {
// printf("lala\\n");
//}
//if (%(z)s) {
// Py_XDECREF(%(z)s);
//}
%(z)s = transposed;
Py_XINCREF(%(z)s);
"""
from gof import modes
modes.make_constructors(globals())
def scalar_switch(normal_f, scalar_f, scalar_f_reverse = None):
def f(x, y):
x, y = wrap_as_tensor(x), wrap_as_tensor(y)
if 0 not in y.broadcastable:
return scalar_f(x, y)
if 0 not in x.broadcastable:
if scalar_f_reverse:
return scalar_f_reverse(y, x)
else:
raise TypeError("You cannot do this operation on a scalar.")
return normal_f(x, y)
return f
# Wrapper to ensure that all inputs to the function impl have the same size (foils numpy's broadcasting)
def assert_same_shapes(x, *rest):
shape = x.shape
for other in rest:
if other.shape != shape:
raise ValueError("The dimensions of the inputs do not match.")
# Wrapper to ensure that the last input to impl is a scalar
def assert_tensor_scalar(x, a):
if numpy.product(a.shape) != 1:
raise ValueError("The second argument must be a scalar.")
class tensor_scalar_op(elemwise):
@classmethod
def variable_names(cls):
return (['x', '_a'], ['z', ])
@classmethod
def loop_variables(cls):
return (['x', ], ['z', ])
def c_init((x, _a), (z, )):
return """
if (PyArray_SIZE(_a) != 1) {
PyErr_SetString(PyExc_ValueError, \"The size of the scalar argument is not 1.\");
}
_a_dtype a = ((_a_dtype*)PyArray_DATA(_a))[0];
"""
def _c_foreach(self):
return "z_i = %s;" % self.c_expr
#########
## Add ##
#########
# Elemwise #
class AddElemwise(Elemwise):
def impl(self, x, y):
assert_same_shapes(x, y)
return x + y
def grad(self, (x, y), gz):
return gz, gz
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i + y_i;"
class AddElemwiseInplace(AddElemwise.inplace_version()):
def impl(self, x, y):
assert_same_shapes(x, y)
x += y
return x
# Scalar #
class AddScalar(TensorScalarOp):
def impl(self, x, a):
assert_tensor_scalar(x, a)
return x + a
def grad(self, (x, a), gz):
return gz, sum(gz)
c_expr = "x_i + a"
AddScalarInplace = add_scalar.inplace_version()
add_scalar_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__iadd__))
class twice(elemwise):
def impl(x):
return 2.0 * x
def grad(x, gz):
return scale(gz, 2.0)
def c_foreach((x_i, ), (z_i, )):
"z_i = x_i + x_i;"
twice_inplace = twice.inplace_version()
## Subtraction ##
class sub_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__sub__)
def grad(x, y, gz):
return gz, -gz
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i - y_i;"
sub_elemwise_inplace = sub_elemwise.inplace_version()
sub_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__isub__))
def sub_scalar_r(x, a):
return add_scalar(x, -a)
def sub_scalar_l(x, a):
return add_scalar(-x, a)
def sub_scalar_r_inplace(x, a):
return add_scalar_inplace(x, -a)
## Element-wise multiplication ##
class mul_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__mul__)
def grad(x, y, gz):
return mul(y, gz), mul(x, gz)
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i * y_i;"
mul_elemwise_inplace = mul_elemwise.inplace_version()
mul_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__imul__))
class scale(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__mul__)
def grad(x, a, gz):
return scale(a, gz), sum(mul_elemwise(x, gz))
c_expr = "x_i * a"
scale_inplace = scale.inplace_version()
scale_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__imul__))
class sqr(elemwise):
def impl(x):
return x * x
def grad(x, gz):
return scale(mul_elemwise(x, gz), 2.0)
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i * x_i;"
isqr = sqr.inplace_version()
isqr.set_impl(lambda x: x.__imul__(x))
class sqrt(elemwise):
impl = numpy.sqrt
def grad(x, gz):
return scale(div(gz, sqrt(x)), 0.5)
def c_foreach((x_i, ), (z_i, )):
return "z_i = pow(x_i, 0.5);"
isqrt = sqrt.inplace_version()
isqrt.set_impl(lambda x: x.__ipow__(0.5))
## Element-wise division ##
class div_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__div__)
def grad(x, y, gz):
return div(gz, y), -div(mul(x, gz), sqr(y))
def c_foreach((x_i, y_i), (z_i, )):
return "z_i = x_i / y_i;"
div_elemwise_inplace = div_elemwise.inplace_version()
div_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__idiv__))
def div_scalar_r(x, a):
return scale(x, inv_elemwise(a))
def div_scalar_l(x, a):
return scale(inv_elemwise(x), a)
def div_scalar_r_inplace(x, a):
return scale_inplace(x, inv_elemwise(a))
## Scaling ##
class scale(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__mul__)
def grad(x, a, gz):
return scale(a, gz), sum(mul_elemwise(x, gz))
c_expr = "x_i * a"
scale_inplace = scale.inplace_version()
scale_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__imul__))
class neg(elemwise):
impl = numpy.ndarray.__neg__
def grad(x, gz):
return -gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = -x_i;"
neg_inplace = neg.inplace_version()
neg_inplace.set_impl(lambda x: x.__imul__(-1))
class inv_elemwise(elemwise):
impl = lambda x: 1 / x
def grad(x, gz):
return -gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = 1 / x_i;"
inv_elemwise_inplace = inv_elemwise.inplace_version()
## Dot product ##
class dot(omega_op):
@staticmethod
def _output_shape(xshape, yshape):
if len(xshape) == 0: # x is a scalar
shape = yshape
else:
if len(yshape) >= 2: #y is a matrix or tensor
assert xshape[-1] == yshape[-2]
shape = tuple(xshape[:-1]+ yshape[:-2]+yshape[-1:])
elif len(yshape)==1: #y is vector
assert xshape[-1] == yshape[-1]
shape = tuple(xshape[:-1])
else: #y is a scalar
shape = xshape
return shape
impl = numpy.dot
def grad(x, y, gz):
return dot(gz, transpose(y)), dot(transpose(x), gz)
def refresh(self, alloc=False):
x,y = self.inputs
shape = self._output_shape(x.shape, y.shape)
dtype = upcast(x.dtype, y.dtype)
if self.out.data is not None \
and self.out.shape == shape \
and self.out.dtype == dtype:
return #everything is ok
if alloc or self.out.data is not None: #data should be allocated
self.out.data = None
self.out.shape = shape
self.out.dtype = dtype
self.out.alloc()
else:
self.out.shape = shape
self.out.dtype = dtype
def c_support_code(self):
return blas.cblas_header_text()
def c_libs(self):
return blas.ldflags()
def c_impl((_x, _y), (_z, )):
return blas.gemm_code('', '1.0', '0.0')
## Transposition ##
class transpose(omega_op):
def view_map(self): return {self.out: [self.inputs[0]]}
impl = numpy.transpose
def grad(x, gz):
return transpose_copy(gz)
def refresh_shape(self):
rval = list(self.inputs[0].shape)
rval.reverse()
return rval
def refresh_dtype(self):
return self.inputs[0].dtype
def c_impl((x, ), (xt, )):
return """
const int l = x->nd;
// The user must ensure that all references to
//xt->data go through xt, or there's going to be trouble..
int refcheck = 0;
if (x == xt)
{
return -1;
}
if (refcheck)
{
int refcnt = PyArray_REFCOUNT(xt);
if ((refcnt > 2) // you might think this should be 1.. but this works
//|| (xt->base != NULL)
|| (xt->weakreflist != NULL))
{
PyErr_SetString(PyExc_ValueError,
"cannot resize an array that has "\\
"been referenced or is referencing\\n"\\
"another array in this way. Use the "\\
"resize function");
return -2;
}
}
if (xt->nd != x->nd)
{
// this technique comes from PyArray_Resize()
npy_intp * dimptr = (npy_intp*)PyDimMem_RENEW(xt->dimensions, 2 * x->nd);
if (!dimptr)
{
PyErr_NoMemory();
return 1;
}
xt->nd = x->nd;
xt->dimensions = dimptr;
xt->strides = dimptr + x->nd;
}
//copy x's dimensions and strides
for (int i = 0; i < l; ++i)
{
xt->dimensions[i] = x->dimensions[l-i-1];
xt->strides[i] = x->strides[l-i-1];
}
// point directly at b's type descriptor
Py_INCREF(x->descr);
Py_DECREF(xt->descr);
xt->descr = x->descr;
// name x as a base of xt, increment its refcount
if ( xt->base != (PyObject*)x)
{
Py_INCREF(x);
if ((xt->base) && (xt->base != Py_None))
{
Py_DECREF(xt->base);
}
xt->base = (PyObject*)x;
}
// mark xt as not owning its data
if (PyArray_CHKFLAGS(xt, NPY_OWNDATA))
{
PyDataMem_FREE(xt->data);
xt->flags &= ~NPY_OWNDATA;
}
xt->data = x->data;
// this function is described in
// ~/zzz.NOBACKUP/pub/src/numpy-1.0.3.1/numpy/core/src/arrayobject.c:1890
PyArray_UpdateFlags(xt, NPY_CONTIGUOUS|NPY_FORTRAN|NPY_ALIGNED|NPY_WRITEABLE);
/*
TODO
What should be done with the weakreflist ?
*/
"""
def transpose_copy(x):
return array_copy(transpose(x))
## Copy ##
class array_copy(elemwise):
impl = numpy.array
grad = lambda x, gz: gz
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i;"
## Power ##
class sqr(elemwise):
def impl(x):
return x * x
def grad(x, gz):
return scale(mul_elemwise(x, gz), 2.0)
def c_foreach((x_i, ), (z_i, )):
return "z_i = x_i * x_i;"
sqr_inplace = sqr.inplace_version()
sqr_inplace.set_impl(lambda x: x.__imul__(x))
class sqrt(elemwise):
impl = numpy.sqrt
def grad(x, gz):
return scale(div(gz, sqrt(x)), 0.5)
def c_foreach((x_i, ), (z_i, )):
return "z_i = pow(x_i, 0.5);"
sqrt_inplace = sqrt.inplace_version()
sqrt_inplace.set_impl(lambda x: x.__ipow__(0.5))
class exp(elemwise):
def impl(x): return numpy.exp(x)
def grad(x, gz): return gz * exp(x)
def c_foreach((x_i, ), (z_i, )): return "z_i = exp(x_i);"
class log(elemwise):
def impl(x): return numpy.log(x)
def grad(x, gz): return gz / x
def c_foreach((x_i, ), (z_i, )): return "z_i = log(x_i);"
class log2(elemwise):
def impl(x): return numpy.log2(x)
def grad(x, gz): return gz / (x * numpy.log(2))
def c_foreach((x_i, ), (z_i, )): return "z_i = log2(x_i);"
class pow_elemwise(elemwise):
impl = assert_same_shapes(numpy.ndarray.__pow__)
def grad(x, s, gz):
raise NotImplemented # no gs
return gz * s * (pow_elemwise(x, s-1.0))
def c_foreach((x_i, s_i), (z_i, )):
return "z_i = pow(x_i, s_i)"
pow_elemwise_inplace = pow_elemwise.inplace_version()
pow_elemwise_inplace.set_impl(assert_same_shapes(numpy.ndarray.__ipow__))
class pow_scalar_l(tensor_scalar_op):
impl = tensor_scalar_impl(lambda x, y: numpy.ndarray.__pow__(y, x))
def grad(x, s, gz):
raise NotImplemented # no gs
return gz * x * (pow_scalar_l(s,x-1.0))
c_expr = "pow(a, x_i)"
class pow_scalar_r(tensor_scalar_op):
impl = tensor_scalar_impl(numpy.ndarray.__pow__)
def grad(x, s, gz):
gx = gz * s * (pow_scalar_r(x,s-1.0))
gs = sum(gz * pow_scalar_r(x,s) * log(x))
return gx, gs
c_expr = "pow(x_i, a)"
pow_scalar_r_inplace = pow_scalar_r.inplace_version()
pow_scalar_r_inplace.set_impl(tensor_scalar_impl(numpy.ndarray.__ipow__))
## Others ##
class minmax(elemwise):
nout = 2
def impl(x):
return x.min, x.max
def specs(x):
return [(numpy.ndarray, x[1], ())] * 2
# def alloc((x, ), (_min, _max)):
# _min.data = numpy.ndarray((), x.dtype)
# _max.data = numpy.ndarray((), x.dtype)
def c_init((x, ), (_min, _max)):
raise NotImplementedError
return """
_x_dtype min = _x[0];
_x_dtype max = _x[0];
"""
def c_foreach((x, ), (_min, _max)):
return """
if (x < min) min = x;
if (x > max) max = x;
"""
def c_finalize((x, ), (_min, _max)):
return """
_min[0] = min;
_max[0] = max;
"""
class fill(elemwise):
impl = lambda model, value: (model * 0) + value
def c_init((model, value), (z, )):
return "value_dtype value0 = ((value_dtype*)PyArray_DATA(value))[0];"
def c_foreach((model_i, value), (z_i, )):
return "z_i = value0;"
fill_inplace = fill.inplace_version()
class sum(elemwise):
impl = numpy.sum
def grad(x, gz):
return fill(x, gz)
def refresh_shape(self):
return ()
def c_init((x, ), (sum, )):
return "sum_dtype* sump = ((sum_dtype*)PyArray_DATA(sum)); sump[0] = 0;"
def c_foreach((x_i, ), (sum, )):
return "sump[0] += x_i;"
class ones_like(elemwise):
impl = numpy.ones_like
def grad(x, gz): return Undefined
class zeros_like(elemwise):
impl = numpy.zeros_like
def grad(x, gz): return Undefined
## Array slicing ##
class get_slice(omega_op):
def view_map(self): return {self.out: [self.inputs[0]]}
def impl(x, item):
rval = x.__getitem__(item)
#print 'get_slice running', rval
return rval
def grad(x, gz): raise NotImplemented
def refresh_shape(self):
x,item = self.inputs
rval = x.data.__getitem__(item.data).shape
#print 'refresh_shape', rval
return rval
def refresh_dtype(self):
return self.inputs[0].data.dtype
add = scalar_switch(add_elemwise, add_scalar, add_scalar)
add_inplace = scalar_switch(add_elemwise_inplace, add_scalar_inplace)
sub = scalar_switch(sub_elemwise, sub_scalar_r, sub_scalar_l)
sub_inplace = scalar_switch(sub_elemwise_inplace, sub_scalar_r_inplace)
mul = scalar_switch(mul_elemwise, scale, scale)
mul_inplace = scalar_switch(mul_elemwise_inplace, scale_inplace)
div = scalar_switch(div_elemwise, div_scalar_r, div_scalar_l)
div_inplace = scalar_switch(div_elemwise_inplace, div_scalar_r_inplace)
pow = scalar_switch(pow_elemwise, pow_scalar_r, pow_scalar_l)
pow_inplace = scalar_switch(pow_elemwise_inplace, pow_scalar_r_inplace)
from scipy.weave import c_spec, standard_array_spec
class omega_type_converter_extension:
# def provides(self):
# """
# Returns a list of (c_type, name, init_code) tuples that represent variables
# the type converter provides to the user's code.
# """
# return []
# def format_provide(self, x):
# return '%s %s = %s;\n' % x
def declaration_code(self, templatize = 0, inline = 0):
tvars = self.template_vars(inline=inline)
code = '%(py_var)s = %(var_lookup)s;\n' % tvars
code += '%(c_type)s %(name)s = %(var_convert)s;\n' % tvars
return code
def struct_init_code(self):
return "Py_INCREF(py_%s);" % self.name
def struct_cleanup_code(self):
return "Py_DECREF(py_%s);" % self.name
def struct_members_code(self):
tvars = self.template_vars()
res = "PyObject* py_%s;\n" % self.name
res += "%(c_type)s %(name)s;\n" % tvars
return res
def struct_import_code(self):
res = "__STRUCT_P->py_%s = py_%s;\n" % (self.name, self.name)
res += "__STRUCT_P->%s = %s;\n" % (self.name, self.name)
return res
def struct_support_code(self):
return ""
def struct_typedefs(self):
return ""
class int_converter(omega_type_converter_extension, c_spec.int_converter):
pass
class float_converter(omega_type_converter_extension, c_spec.float_converter):
pass
class complex_converter(omega_type_converter_extension, c_spec.complex_converter):
pass
class unicode_converter(omega_type_converter_extension, c_spec.unicode_converter):
pass
# def provides(self):
# tvars = self.template_vars()
# return omega_type_converter_extension.provides() + [('int', 'N%(name)s' % tvars, 'PyUnicode_GET_SIZE(%(py_var)s)' % tvars)]
class string_converter(omega_type_converter_extension, c_spec.string_converter):
pass
class list_converter(omega_type_converter_extension, c_spec.list_converter):
pass
class dict_converter(omega_type_converter_extension, c_spec.dict_converter):
pass
class tuple_converter(omega_type_converter_extension, c_spec.tuple_converter):
pass
class file_converter(omega_type_converter_extension, c_spec.file_converter):
pass
class instance_converter(omega_type_converter_extension, c_spec.instance_converter):
pass
class array_converter(omega_type_converter_extension, standard_array_spec.array_converter):
# def provides(self):
# tvars = self.template_vars()
# ret = []
# ret.append((tvars['c_type'], tvars['array_name'], tvars['var_convert']))
# ret.append(('npy_intp*', 'N%(name)s' % tvars, '%(array_name)s->dimensions' % tvars))
# ret.append(('npy_intp*', 'S%(name)s' % tvars, '%(array_name)s->strides' % tvars))
# ret.append(('int', 'D%(name)s' % tvars, '%(array_name)s->nd' % tvars))
# ret.append(('%(num_type)s*' % tvars, '%(name)s' % tvars, '(%(num_type)s*) %(array_name)s->data' % tvars))
# return ret
# def declaration_code(self, templatize = 0, inline = 0):
# tvars = self.template_vars(inline=inline)
# tvars['cap_name'] = self.name.upper()
# prov = self.provides()
# code = '%(py_var)s = %(var_lookup)s;\n' % tvars
# code += "\n".join(self.format_provide(export) for export in prov[:1])
# code += '\nconversion_numpy_check_type(%(array_name)s,%(num_typecode)s,"%(name)s");\n' % tvars
# code += "\n".join(self.format_provide(export) for export in prov[1:])
# return code
# def struct_support_code(self, templatize = 0, inline = 0):
# tvars = self.template_vars(inline=inline)
# cap_name = self.name.upper()
# tvars['cap_name'] = cap_name
# code = 'inline %(num_type)s& %(cap_name)s1(int i) { return (*((%(num_type)s*)(%(array_name)s->data + (i)*S%(name)s[0])));}\n' \
# 'inline %(num_type)s& %(cap_name)s2(int i, int j) { return (*((%(num_type)s*)(%(array_name)s->data + (i)*S%(name)s[0] + (j)*S%(name)s[1])));}\n' \
# 'inline %(num_type)s& %(cap_name)s3(int i, int j, int k) { return (*((%(num_type)s*)(%(array_name)s->data + (i)*S%(name)s[0] + (j)*S%(name)s[1] + (k)*S%(name)s[2])));}\n' \
# 'inline %(num_type)s& %(cap_name)s4(int i, int j, int k, int l) { return (*((%(num_type)s*)(%(array_name)s->data + (i)*S%(name)s[0] + (j)*S%(name)s[1] + (k)*S%(name)s[2] + (l)*S%(name)s[3])));}\n'
# return code % tvars
def struct_typedefs(self):
tvars = self.template_vars()
return "typedef %(num_type)s %(name)s_dtype;\n" % tvars
# return "\n".join(["typedef %s %s_type;" % (c_type, name)])
# def struct_template_types(self):
# tvars = self.template_vars()
# return [("typename %s_type" % name, c_type) for c_type, name, init in self.provides()] + [("typename %s_dtype" % self.name, tvars['num_type'])]
default = [array_converter(),
int_converter(),
float_converter(),
complex_converter(),
unicode_converter(),
string_converter(),
list_converter(),
dict_converter(),
tuple_converter(),
file_converter(),
instance_converter()]
from core import Numpy2, omega_op
def input(x):
#static member initialization
if not hasattr(input, 'float_dtype'):
input.float_dtype = 'float64'
input.int_dtype = 'int64'
input.NN = Numpy2
if isinstance(x, numpy.ndarray):
#return NumpyR(x)
return input.NN(data=x)
elif isinstance(x, int):
z = numpy.zeros((), dtype = input.int_dtype)
z += x
return input.NN(data=z)
elif isinstance(x, float):
z = numpy.zeros((), dtype = input.float_dtype)
z += x
return input.NN(data=z)
elif is_result(x):
raise TypeError("%s is already a result." % x)
else:
return ResultBase(data=x)
def wrap(x):
if isinstance(x, Numpy2):
return x
#elif isinstance(x, NumpyR):
#return x
elif is_result(x):
return x
elif isinstance(x, omega_op):
return x.out
else:
return literal(x)
def literal(x):
"""Return a ResultValue instance wrapping a literal."""
def _hashable(x):
try:
x in {}
return True
except TypeError: # x is unhashable
return False
#static member initialization
if not hasattr(literal, 'hdb'):
literal.hdb = {}
literal.udb = {}
if _hashable(x):
db = literal.hdb
key = (type(x),x)
else:
db = literal.udb
key = (id(x),)
if key in db:
return db[key]
else:
rval = input(x)
rval.constant = True
db[key] = rval
return rval
Markdown 格式
0%
您添加了 0 到此讨论。请谨慎行事。
请先完成此评论的编辑!
注册 或者 后发表评论