<jupyter_start><jupyter_text>Find initial guess<jupyter_code># Parameters of our force solving method
contactPadding = 15
g2MaskPadding = 3
contactMaskRadius = 30
#pxPerMeter = 10000
brightfield = False

forceGuessArr, alphaGuessArr, betaGuessArr = initialForceSolve(photoelasticChannel,
                                                centers, radii, fSigma, pxPerMeter,
                                                contactPadding, g2MaskPadding,
                                                contactMaskRadius=contactMaskRadius,
                                                boundaryMask=maskArr, ignoreBoundary=False)

#forceGuessArr /= np.max([np.max(f) for f in forceGuessArr])<jupyter_output><empty_output><jupyter_text>Show initial guess<jupyter_code># Draw on all of the found contacts
fig, ax = plt.subplots()

ax.imshow(maskArr*255)

for i in range(len(centers)):
    # particle itself
    c = plt.Circle(centers[i][::-1], radii[i], label='Detected particle', color='orange')
    ax.add_artist(c)
    
    # Now add contacts
    for j in range(len(betaGuessArr[i])):
        contactPoint = centers[i] + radii[i] * np.array([np.cos(betaGuessArr[i][j]), np.sin(betaGuessArr[i][j])])
        #linePoints = np.array([centers[i], radii[i] * np.array([np.cos(betaGuessArr[i][j]), np.sin(betaGuessArr[i][j])])])
        #plt.plot(linePoints[0], linePoints[1], color='blue')
        cc = plt.Circle(contactPoint[::-1], 20, color='red')
        ax.add_artist(cc)
        
plt.show()
fig, ax = plt.subplots(1, 2, figsize=(8,4))

ax[0].imshow(photoelasticChannel * maskArr[:,:,0])
ax[0].set_title('True photoelastic response\n(unknown forces)')

estimatedPhotoelasticChannel = np.zeros_like(photoelasticChannel, dtype=np.float64)

for i in range(len(centers)):
    estimatedPhotoelasticChannel += genSyntheticResponse(np.array(forceGuessArr[i]),
                                                         np.array(alphaGuessArr[i]),
                                                         np.array(betaGuessArr[i]),
                                                         fSigma, radii[i],
                                                         pxPerMeter, brightfield, imageSize=photoelasticChannel.shape,
                                                         center=centers[i])


ax[1].imshow(estimatedPhotoelasticChannel)
ax[1].set_title('Initial Guess for Optimizer\n(known forces)')

fig.tight_layout()
plt.show()

print(f'Absolute error (entire system): {np.sum(estimatedPhotoelasticChannel - photoelasticChannel)}')<jupyter_output><empty_output><jupyter_text>Optimize forces<jupyter_code># Now optimize each particle's forces
optimizedForceArr = []
optimizedBetaArr = []
optimizedAlphaArr = []

for i in range(len(centers)):
    optForceArr, optBetaArr, optAlphaArr, res = forceOptimize(forceGuessArr[i], betaGuessArr[i], alphaGuessArr[i], radii[i], centers[i], photoelasticChannel,
                                                              fSigma, pxPerMeter, brightfield, maxEvals=[250, 150], method='nelder',
                                                              parametersToFit=[['f', 'a'], ['f', 'b']], betaTol=.4, forceTol=5., allowAddForces=True)
    optimizedForceArr.append(optForceArr)
    optimizedBetaArr.append(optBetaArr)
    optimizedAlphaArr.append(optAlphaArr)
fig, ax = plt.subplots(1, 2, figsize=(8,4))

ax[0].imshow(photoelasticChannel * maskArr[:,:,0])
ax[0].set_title('True photoelastic response\n(unknown forces)')
ax[0].set_xticks([])
ax[0].set_yticks([])

optimizedPhotoelasticChannel = np.zeros(photoelasticChannel.shape)

for i in range(len(centers)):
    optimizedPhotoelasticChannel += genSyntheticResponse(np.array(optimizedForceArr[i]),
                                                         np.array(optimizedAlphaArr[i]),
                                                         np.array(optimizedBetaArr[i]),
                                                         fSigma, radii[i], pxPerMeter, brightfield,
                                                         imageSize=photoelasticChannel.shape, center=centers[i])


ax[1].imshow(optimizedPhotoelasticChannel)
ax[1].set_title('Optimized photoelastic response\n(known forces)')
ax[1].set_xticks([])
ax[1].set_yticks([])

fig.tight_layout()
#plt.savefig('../images/example_uses/space_jam_test.png')
plt.show()

print(f'Absolute error (entire system): {np.sum(optimizedPhotoelasticChannel - photoelasticChannel)}')
optimizedForceArr
# To try and identify missing forces, take the average g2 in a ring
# around the edge of a particle
testParticle = 0
contactMaskRadius = 30.
steps = 30
newBetaTolerance = .3

avgG2Arr = np.zeros(steps)
betaArr = np.linspace(-np.pi, np.pi, steps)

gSqr = gSquared(photoelasticChannel)
particleMask = circularMask(photoelasticChannel.shape, centers[testParticle], radii[testParticle])[:,:,0]

for i in range(steps):
    contactPoint = centers[testParticle] + radii[testParticle] * np.array([np.cos(betaArr[i]), np.sin(betaArr[i])])
    
    contactMask = circularMask(photoelasticChannel.shape, contactPoint, contactMaskRadius)[:,:,0]
    contactMask = (contactMask + particleMask) == 2
    
    avgG2Arr[i] = np.sum(contactMask * gSqr) / np.sum(contactMask)
    
peakIndArr = find_peaks(avgG2Arr, height=.0003)[0]
peakIndArr = np.sort(peakIndArr)

# Make sure that there aren't any artifacts of periodicity
if np.arange(3).any() in peakIndArr and np.arange(len(avgG2Arr)-3, len(avgG2Arr)).any() in peakIndArr:
    print('overlap')
    # Remove last entry
    peakIndArr = peakIndArr[:-1]
    
peakBetas = betaArr[peakIndArr]

plt.plot(betaArr, avgG2Arr)

for i in range(len(peakBetas)):
    plt.axvline(peakBetas[i], linestyle='--')

plt.show()
# Now determine which of the found peaks is the most different
# from the ones already on the particle
oldBetaArr = betaGuessArr[testParticle]

differenceArr = np.abs(np.subtract.outer(peakBetas, oldBetaArr))

# Check to see if there is a new peak that doesn't have
# a previous force close to it
for i in range(len(peakBetas)):
    if differenceArr[i].all() > newBetaTolerance:
        # Add the new force
# How do we identify when we need to add a force?
testParticle = 3
result = forceOptimize(forceGuessArr[testParticle], betaGuessArr[testParticle], alphaGuessArr[testParticle], radii[testParticle], centers[testParticle], photoelasticChannel,
                                        fSigma, pxPerMeter, brightfield, maxEvals=150, method='nelder',
                                        parametersToFit=['f', 'b'], betaTol=.4, returnOptResult=True)
print(fit_report(result))
result.chisqr<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1202416918429003, 'fraction_numerical': 0.012386706948640483, 'mean_word_length': 2.9551971326164876, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '636633', 'license': 'MIT', 'n_tokens_mistral': 2101, 'n_tokens_neox': 1995, 'n_words': 400, 'path': 'notebooks/force_adding_test.ipynb', 'repo_name': 'Jfeatherstone/pepe', 'size': '6620'}