Given a function \(\Phi(x,n) = \left|\phi(x+n)-\phi(x)-\phi(n)\right| \), the algorithm below will find \( \phi(x) \) using \( \Phi(x,1) \) and \( \Phi(x,2) \), \( \phi(1) \) and \( \phi(2) \) and \( \phi(3) \), and the anti-symmetry of \( \phi(x) \) (which also gives us \( \phi(0) = 0 \)). This can be extended so that the algorithm uses \( \Phi(x,n) \) for more values of \( n \) to prevent mistakes in the choice of path for certain critical forms of \( \phi(x) \) but those cases seldomly appear.
I was not able to figure out a method to reconstruct \( \phi(x) \) without knowing the first three non-trivial values; the anti-symmetry defeats you in this case, I think. Probably these values can be constrained and a decent enough optimization made of the entire function, but this is not ideal.
import numpy as np
import matplotlib.pyplot as plt
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return array[idx]
def PhiSolve(size = 10):
tmp = np.random.random(size=size)*20 - 10
phi = np.concatenate((-tmp, [0], np.flip(tmp))).astype(int)
PHI_1 = np.abs(np.roll(phi, -1) - phi - phi[len(phi)//2+1])
PHI_2 = np.abs(np.roll(phi, -2) - phi - phi[len(phi)//2+2])
# Algorithm
solved = np.zeros_like(phi)
# Even if we don't know phi(2) and phi(3), calculating for those will give two options each for a total of four possible branches
# Still need to know phi(0) (trivial) and phi(1) to get the four choices. There are an infinite number of naive choices for phi(1), however...
solved[len(phi)//2+1] = phi[len(phi)//2+1]
solved[len(phi)//2+2] = phi[len(phi)//2+2]
solved[len(phi)//2+3] = phi[len(phi)//2+3]
for n in range(len(solved)//2+3,len(solved)-1,1):
x1 = PHI_1[n] + solved[n] + solved[len(solved)//2+1]
x2 = -PHI_1[n] + solved[n] + solved[len(solved)//2+1]
x3 = PHI_2[n-1] + solved[n-1] + solved[len(solved)//2+2]
x4 = -PHI_2[n-1] + solved[n-1] + solved[len(solved)//2+2]
sorted = np.sort([x1,x2,x3,x4])
diff = np.abs([sorted[0]-sorted[1], sorted[1]-sorted[2],sorted[2]-sorted[3]])
idx = np.argmin(diff)
solved[n+1] = np.mean([sorted[idx],sorted[idx+1]])
solved[:len(solved)//2] = np.flip(-solved[len(solved)//2+1:])
plt.plot(np.linspace(-size,size,len(phi)), solved)
plt.plot(np.linspace(-size, size, len(phi)), phi, 'o--')
plt.show()
plt.close()
print(phi)
print(PHI_1)
print(PHI_2)
if __name__ == '__main__':
PhiSolve(size=10)
