stackexchange

import numpy as np
import scipy.optimize as so

def synth_I_field_compact(part_params, nd, sigma, x, y, z=np.zeros(0), tol=1e-4):
    nparts = int(np.size(part_params) / (nd + 1))
    xpart = part_params[0:nparts] # must be column vector from reshape
    ypart = part_params[nparts:2 * nparts]
    zpart = part_params[nparts:3 * nparts]
    s = part_params[nd*nparts:]
    nx = np.size(x)
    ny = np.size(y)
    nz = np.size(z)
    I = np.zeros((ny, nx, nz))
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    dz = z[1] - z[0]
    di = np.ceil(max_dist/dx)

    for i in range(0,nparts):
        ix = round((xpart[i] - x[0]) / dx)
        iy = round((ypart[i] - y[0]) / dy)
        iz = round((ypart[i] - y[0]) / dz)
        iix = np.arange(max(0, ix - di), min(nx, ix + di), 1, dtype=int)
        iiy = np.arange(max(0, iy - di), min(ny, iy + di), 1, dtype=int)
        iiz = np.arange(max(0, iz - di), min(nz, iz + di), 1, dtype=int)
        ddx = dx * iix - xpart[i]
        ddy = dy * iiy - ypart[i]
        ddz = dz * iiz - zpart[i]
        gx = np.exp(-1 / (2 * sigma ** 2) * ddx ** 2)
        gy = np.exp(-1 / (2 * sigma ** 2) * ddy ** 2)
        gz = np.exp(-1 / (2 * sigma ** 2) * ddz ** 2)
        gx = gx[np.newaxis,:, np.newaxis]
        gy = gy[:,np.newaxis, np.newaxis]
        gz = gz[np.newaxis, np.newaxis, :]
        I[np.ix_(iiy, iix, iiz)] = I[np.ix_(iiy, iix, iiz)] + s[i] * gy*gx*gz
    return I

nd = 3
x = np.linspace(0, 1, 29) # spatial domain
y = np.linspace(0, 1, 31)
z = np.linspace(0, 1, 30)

xp = np.asarray([0.5]) # particle position
yp = np.asarray([0.5])
zp = np.asarray([0.5])

nx = np.size(x)
ny = np.size(y)
nz = np.size(z)

I = np.zeros((ny, nx, nz)) # intensity field
dx = x[1] - x[0]
dy = y[1] - y[0]
dz = z[1] - z[0]

tolerance = 1e-4 # minimum value of Gaussian for compact representation
sigma = x[1]-x[0] # Gaussian width
max_dist = sigma*(-2*np.log(tolerance))
di = np.ceil(max_dist/dx) # maximum distance in compact representation, in index format

# Create intensity field/true Gaussian
# this exists separately as its own function synth_I() where [0] is instead for each particle [i]
ix = round((xp[0] - x[0]) / dx) # index of particle center
iy = round((yp[0] - y[0]) / dy)
iz = round((yp[0] - y[0]) / dz)
iix = np.arange(max(0, ix - di), min(nx, ix + di), 1, dtype=int) # grid points with nonzero intensity values
iiy = np.arange(max(0, iy - di), min(ny, iy + di), 1, dtype=int)
iiz = np.arange(max(0, iz - di), min(nz, iz + di), 1, dtype=int)
ddx = dx * iix - xp[0] # distance between particle center and grid point
ddy = dy * iiy - yp[0]
ddz = dz * iiz - zp[0]
gx = np.exp(-1 / (2 * sigma ** 2) * ddx ** 2) # 1D Gaussian
gy = np.exp(-1 / (2 * sigma ** 2) * ddy ** 2)
gz = np.exp(-1 / (2 * sigma ** 2) * ddz ** 2)
gx = gx[np.newaxis,:, np.newaxis]
gy = gy[:,np.newaxis, np.newaxis]
gz = gz[np.newaxis, np.newaxis, :]
I[np.ix_(iiy, iix, iiz)] = I[np.ix_(iiy, iix, iiz)] + gy*gx*gz

# Idea: fit Gaussians centered at grid point to unknown distribution, by adjusting intensity of each
Imin = 0.5
ny, nx, nz = I.shape
X, Y, Z = np.meshgrid(x, y, z)
xpart = X[I > Imin] # estimated particle positions based on minimum intensity threshold
ypart = Y[I > Imin]
zpart = Z[I > Imin]
s_init = I[I> Imin] # initial intensity guess to be used in optimization
Np = np.size(xpart) # number of particles to fit to
# simulate intensity as sum of gaussians.

# fit using conjugate gradient algorithm as opposed to least squares, for scalability in real problem
# Objective function: 0.5||I - \sum s[i] * G(x;xpart[i])||^2
# Gradient: g[j] =  - \sum { I - \sum s[i] * G(x;xpart[i]) * G(x;xpart[j]) }

def diff_and_grad(s):  # ping test this:
    part_params = np.concatenate((xpart, ypart, zpart, s)) # for own code
    # create an intensity field as a combination of Gaussians as above
    synthI = synth_I_field_compact(part_params, nd, sigma, x, y, z)
    Idiff = I - synthI # difference in measurements
    f = 0.5 * np.sum(np.power(Idiff, 2)) # objective function
    g = np.zeros(Np) # gradient
    for i in range(0, Np):
        ix = round((xpart[i] - x[0]) / dx)
        iy = round((ypart[i] - y[0]) / dy)
        iz = round((zpart[i] - z[0]) / dz)
        iix = np.arange(max(0, ix - di), min(nx, ix + di), 1, dtype=int)
        iiy = np.arange(max(0, iy - di), min(ny, iy + di), 1, dtype=int)
        iiz = np.arange(max(0, iz - di), min(nz, iz + di), 1, dtype=int)
        ddx = dx * iix - xpart[i]
        ddy = dy * iiy - ypart[i]
        ddz = dz * iiz - zpart[i]
        gx = np.exp(-1 / (2 * sigma ** 2) * ddx ** 2)
        gy = np.exp(-1 / (2 * sigma ** 2) * ddy ** 2)
        gz = np.exp(-1 / (2 * sigma ** 2) * ddz ** 2)
        Id = Idiff[np.ix_(iiy, iix, iiz)]
        g[i] = -Id.dot(gz).dot(gx).dot(gy) # gradient is -product of local intensity difference with gaussian centered at grid point
    return f, g

# only objective function:
def diff(s):
    part_params = np.concatenate((xpart,ypart,zpart,s))
    synthI = synth_I_field_compact(part_params, nd, sigma, x,y,z)
    Idiff = I - synthI
    f = 0.5*np.sum(np.power(Idiff,2))
    return f

# differing gradients:
f_init, g_ana = diff_and_grad(s_init)
g_fd = np.zeros(Np)
delta = 1e-5
for j in range(0, Np):
    ej = np.zeros(Np)
    ej[j] = 1
    fp = diff(s_init+delta*ej)
    fm = diff(s_init-delta*ej)
    g_fd[j] = 1 / (2 * delta) * (fp - fm)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(g_ana)
plt.hold(True)
plt.plot(g_fd)
plt.legend(["analytical", "finite difference"])
plt.hold(False)
plt.show()

result = so.minimize(diff_and_grad, s_init, method='cg', jac=True)
result2 = so.minimize(diff, s_init, method='cg')