import numpy as np

def MLEM(G, G1, sino, imold, nx, na, nb, R=None, beta=1e-3):
    sino1=G.dot(imold.ravel()).reshape(na,nb)
    sino2=sino.copy()
    sino2 /= sino1
    #G1=np.transpose(G)
    #im1=G1.dot(sino.ravel()).reshape(nx,nx)
    #add regularization
    if R==None:
        im1=G1.dot(sino2.ravel()).reshape(nx,nx)
    else:
        v0=G1.dot(sino2.ravel())
        v1=R.dot(imold.ravel())
        v1=[1+beta*x for x in v1]
        v1=v1.reshape(nx,nx)
        v2=[x/y for x,y in join(v0,v1)]
        im1=v2.reshape(nx,nx)

    im1 *= imold
    return im1