<jupyter_start><jupyter_text>Testing the Mesh Generation Make a mesh with pygmsh, with some dummy values to play with:<jupyter_code>with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    ellipsoid = geom.add_ball([0.0, 0.0, 0.0], 1)
    InputMesh = geom.generate_mesh()
interiorbonds,edgebonds,boundarytris, bidxTotidx, tetras= MakeMeshData3D(InputMesh)
bonds=np.concatenate((interiorbonds,edgebonds))
orientedboundarytris=OrientTriangles(InputMesh.points,boundarytris,np.array([0,0,0]))
cells=[ ("line", bonds ), ("triangle",boundarytris ), ("tetra",tetras)]
isbond=  np.ones(len(bonds))
isedgebond= np.concatenate( ( np.zeros(len(interiorbonds)),np.ones(len(edgebonds)) ) )
CellDataDict={'isedgebond':[isedgebond,np.zeros(len(boundarytris)),np.zeros(len(tetras))]
              ,'isbond':[isbond,np.zeros(len(boundarytris)),np.zeros(len(tetras))]}

OutputMesh=meshio.Mesh(InputMesh.points, cells, {},CellDataDict)
OutputMesh.write(DataFolder+"InitialMesh.vtk",binary=True)<jupyter_output><empty_output><jupyter_text>Testing the boundary triangle finding when pygmsh generates a sphere, it gives the tetrahedrons and boundary triangles. We can use this to check our bondary finding is working. First, make the ball. Now, compare the lists of boundary triangles that we find to those that pygmsh finds. We need to sort the pygmsh ones, as the vertices dont always appear in ascending order<jupyter_code>np.array_equal(boundarytris,np.sort(InputMesh.cells[1].data,axis=1))
boundarytris.shape
InputMesh.cells[1].data.shape<jupyter_output><empty_output><jupyter_text>This seems a good verification that our boundary list is correct  Check total surface area For a sphere<jupyter_code>vTotalArea3D(InputMesh.points,boundarytris)
4*np.pi<jupyter_output><empty_output><jupyter_text>Check total volume Ive written two functions which do this, lets check both<jupyter_code>Volume3D(InputMesh.points,orientedboundarytris,bidxTotidx).sum()
Volume3D_tetras(InputMesh.points,tetras).sum()
(4/3)*np.pi<jupyter_output><empty_output><jupyter_text>Physical Tests  Checking the Bending Modulus Energy  As implemented, the bending modulus approximates the continuum limit $F= \frac{\kappa_c}{2}\int dA(C_1+C_2-C_0)^2 + k_g \int dA C_1C_2$ for a closed surface, where $C_1$ etc. are the principal curvatures. According to Boal and Rao 1992, the energy of a sphere  without a spontaneous curvature is $\frac{4\pi k_{rig}}{\sqrt{3}}$, where $k_{rig}$ is the microscopic modulus.  This seems to be true for a triangulation by equilateral triangles. Lets check this:<jupyter_code>theta_0=0
kbend=1
energies= BendingEnergy(InputMesh.points,orientedboundarytris,bidxTotidx,kbend)
energies.sum()
(4/np.sqrt(3))*np.pi<jupyter_output><empty_output><jupyter_text>Checking the Spring Energy A basic test: lets just make 2 springs (to check the vectorization), and confirm their behaviour. We are supposed to be implementing:$V(r,r_0) = k_{\mathrm{neo}}\left( \frac{1-\alpha}{2}\left(\frac{1}{\lambda}+2\lambda^2\right)+\frac{\alpha}{2}\left(\frac{1}{\lambda^2}+2\lambda\right) \right)$, where $V=r/r_0$, and $k_{neo}=\frac{r_0^2 k_{hook}}{3}$. Some tests:$V(1)= \frac{1}{2} k_{hook} r_0^2=\frac{3}{2} k_{neo}$, indepdent of the mat non.<jupyter_code>MatNon=1
khook=1
#rest lengths
r0_ij=np.array([1,2])
SpringEnergy=NeoHookean3D(r0_ij,r0_ij,khook,MatNon)
print(SpringEnergy) 
np.log(0.5)<jupyter_output><empty_output><jupyter_text>Lets do the same with our shifted energy:<jupyter_code>MatNon=1
khook=1
#rest lengths
r0_ij=np.array([1,2])
SpringEnergy=NeoHookeanShifted(r0_ij,r0_ij,khook,MatNon)
print(SpringEnergy)<jupyter_output><empty_output><jupyter_text>Plotting both energies:<jupyter_code>lam=np.arange(0.1, 10, 0.01);
MatNon=0.7
khook=2

Energy=NeoHookean3D(lam,1,khook,MatNon)
plt.plot(lam,Energy)

kneo_ij = (1**2)*khook/3 
Energy=NeoHookeanShifted(lam,1,khook,MatNon)+1.5*kneo_ij
plt.plot(lam,Energy)<jupyter_output><empty_output><jupyter_text>Plotting the energy: Pure Neohookean on a loglogplot. Expectation: Minimum at $(0,log(0.5)$=$(0,-0.693)$ Asmptotes to grad -1 and grad 2 in either limit.Pure MR: same minimum, but the opposite grad behaviours<jupyter_code>lam=np.arange(0.1, 10, 0.01);
MatNon=0
khook=1
Energy=NeoHookean3D(lam,1,khook,MatNon)
plt.plot(np.log(lam),np.log(Energy))

MatNon=1
khook=1
Energy=NeoHookean3D(lam,1,khook,MatNon)
plt.plot(np.log(lam),np.log(Energy))
lam=10**10
MatNon=0
khook=1
print(( np.log(NeoHookean3D(lam,1,khook,MatNon))-np.log(NeoHookean3D(1,1,khook,MatNon)) )/(np.log(lam)))
lam=10**(-10)
MatNon=0
khook=1
print(( np.log(NeoHookean3D(lam,1,khook,MatNon))-np.log(NeoHookean3D(1,1,khook,MatNon)) )/(np.log(lam)))

lam=10**10
MatNon=1
khook=1
print(( np.log(NeoHookean3D(lam,1,khook,MatNon))-np.log(NeoHookean3D(1,1,khook,MatNon)) )/(np.log(lam)))
lam=10**(-10)
MatNon=1
khook=1
print(( np.log(NeoHookean3D(lam,1,khook,MatNon))-np.log(NeoHookean3D(1,1,khook,MatNon)) )/(np.log(lam)))

matplotlib.rcParams.update({'font.size': 20})
lam=np.arange(0.5,2 , 0.01);
MatNon=0
khook=1
Energy=NeoHookean3D(lam,1,khook,MatNon)
plt.plot((lam),(Energy),linewidth=4.0)

MatNon=1
khook=1
Energy=NeoHookean3D(lam,1,khook,MatNon)
plt.plot((lam),(Energy),linewidth=4.0)

plt.xlabel('$\lambda$')
plt.ylabel('F')
plt.legend(['NH: $\\alpha/\mu=0$','MR: $\\alpha/\\mu=1$'],prop={'size':20})
plt.savefig("Energies.png", bbox_inches='tight',dpi=400)<jupyter_output><empty_output><jupyter_text>Checking that the Numba versions of functions match the regular ones Lets generate a problem similar to our own:<jupyter_code># Target mesh size:
target_a = 0.2
# continuum bending modulus:
kc=0.5
# continuum shear modulus:
mu=1
# Energetic penalty for volume change
B=100000
# The Material Nonlinearity parameter, between 0 and 1
MatNon=0.99
# the spring prestress values 
g0=1

# The microscopic values
kbend=kc/target_a
khook = mu
theta0=0.2
with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = target_a
    ellipsoid = geom.add_ball([0.0, 0.0, 0.0], 1)
    #ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [0.95, 0.95, 1.05])
    InputMesh = geom.generate_mesh()
interiorbonds,edgebonds,boundarytris, bidxTotidx, tetras= MakeMeshData3D(InputMesh)
bonds=np.concatenate((interiorbonds,edgebonds))
orientedboundarytris=OrientTriangles(InputMesh.points,boundarytris,np.array([0,0,0]))
boundarytris=orientedboundarytris
# make the preferred rest lengths of the interior springs
interiorpairs=InputMesh.points[interiorbonds]
interiorvecs = np.subtract(interiorpairs[:,0,:],interiorpairs[:,1,:])
InteriorBondRestLengths=np.linalg.norm(interiorvecs,axis=1)

# make the preferred rest lengths of the edge springs. Initially have the at g0=1, but then
#update them in the loop
edgepairs=InputMesh.points[edgebonds]
edgevecs = np.subtract(edgepairs[:,0,:],edgepairs[:,1,:])
InitialEdgeBondRestLengths=np.linalg.norm(edgevecs,axis=1)

# The volume constraint is simply that the target volume should be the initial volume
TargetVolumes=Volume3D_tetras(InputMesh.points,tetras)

P =InputMesh.points

 # the important bit! Giving it the prestress
EdgeBondRestLengths= g0*InitialEdgeBondRestLengths
r0_ij=np.concatenate((InteriorBondRestLengths,EdgeBondRestLengths))<jupyter_output><empty_output><jupyter_text>To test numerical equality, we can use numpy's testing module:https://numpy.org/doc/stable/reference/generated/numpy.testing.assert_allclose.htmlnumpy.testing.assert_allclose<jupyter_code>x=Volume3D_tetras(P,tetras)
Numbax=NumbaVolume3D_tetras_2(P,tetras)
np.testing.assert_allclose(x,Numbax)
x=BendingEnergy(P,orientedboundarytris,bidxTotidx,kbend)
Numbax=NumbaBendingEnergy_2(P,orientedboundarytris,bidxTotidx,kbend)
np.testing.assert_allclose(x,Numbax)
x=energy3D(P,bonds,orientedboundarytris,bidxTotidx,tetras,r0_ij,khook,kbend,theta0,B,MatNon,TargetVolumes)
Numbax=Numbaenergy3D(P,bonds,orientedboundarytris,bidxTotidx,tetras,r0_ij,khook,kbend,theta0,B,MatNon,TargetVolumes)
np.testing.assert_allclose(x,Numbax)<jupyter_output><empty_output><jupyter_text>Checking Timings  volume<jupyter_code>start = timer()
for i in range(0,5000):
    x=Volume3D_tetras(P,tetras)
end = timer()
print(end-start)
# directly sum the triple product over all tetrahedra
@jit(nopython=True)
def NumbaVolume3D_tetras_2(P,tetras):
    
    Tot=np.zeros(len(tetras))
    for i in range(len(tetras)):
        
        P0= P[tetras[i,0]]
        P1= P[tetras[i,1]] 
        P2= P[tetras[i,2]] 
        P3= P[tetras[i,3]] 
              
        t0x=P1[0]-P0[0]
        t0y=P1[1]-P0[1]
        t0z=P1[2]-P0[2]
        
        t1x=P2[0]-P0[0]
        t1y=P2[1]-P0[1]
        t1z=P2[2]-P0[2]
        
        t2x=P3[0]-P0[0]
        t2y=P3[1]-P0[1]
        t2z=P3[2]-P0[2]
        
        
        t0ct1x = t0y*t1z- t0z*t1y
        t0ct1y = t0z*t1x- t0x*t1z
        t0ct1z = t0x*t1y- t0y*t1x
        
        t2dott0ct1=t2x*t0ct1x+t2y*t0ct1y+t2z*t0ct1z
        
        Tot[i]=np.abs(t2dott0ct1/6)
      
    return Tot
x=NumbaVolume3D_tetras_2(P,tetras)
start = timer()
for i in range(0,5000):
    x=NumbaVolume3D_tetras_2(P,tetras)
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>Bending<jupyter_code>start = timer()
for i in range(0,5000):
    x=BendingEnergy(P,orientedboundarytris,bidxTotidx,kbend)
end = timer()
print(end-start)
x=NumbaBendingEnergy_2(P,orientedboundarytris,bidxTotidx,kbend)
start = timer()
for i in range(0,5000):
    x=NumbaBendingEnergy_2(P,orientedboundarytris,bidxTotidx,kbend)
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>Spring<jupyter_code>start = timer()
for i in range(0,5000):
    x=NeoHookean3D(r0_ij,r0_ij,khook,MatNon).sum()   
end = timer()
print(end-start)
x=NumbaNeoHookean3D(r0_ij,r0_ij,khook,MatNon).sum()  
start = timer()
for i in range(0,5000):
       x=NumbaNeoHookean3D(r0_ij,r0_ij,khook,MatNon).sum()  
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>Making Spring rests<jupyter_code>start = timer()
for i in range(0,5000):
    # We convert it to a matrix here.
    P_ij = P.reshape((-1, 3))
    # from the bond list, work out what the current bond lengths are:
    AB=P_ij[bonds]
    t1 = np.subtract(AB[:,0,:],AB[:,1,:])
    r_ij=np.linalg.norm(t1,axis=1)
end = timer()
print(end-start)
    
# We convert it to a matrix here.
P_ij = P.reshape((-1, 3))
r_ij=NumbaMakeBondLengths(P_ij,bonds)
start = timer()
for i in range(0,5000):
    # We convert it to a matrix here.
    P_ij = P.reshape((-1, 3))
    r_ij=NumbaMakeBondLengths(P_ij,bonds)
end = timer()
print(end-start)
start = timer()
for i in range(0,5000):
    # We convert it to a matrix here.
    P_ij = P.reshape((-1, 3))
    # from the bond list, work out what the current bond lengths are:
    AB=P_ij[bonds]
    t1 = np.subtract(AB[:,0,:],AB[:,1,:])
    r_ij=np.linalg.norm(t1,axis=1)
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>Totals<jupyter_code>start = timer()
for i in range(0,5000):
    x=energy3D(P,bonds,orientedboundarytris,bidxTotidx,tetras,r0_ij,khook,kbend,theta0,B,MatNon,TargetVolumes)
end = timer()
print(end-start)
Numbax=Numbaenergy3D(P,bonds,orientedboundarytris,bidxTotidx,tetras,r0_ij,khook,kbend,theta0,B,MatNon,TargetVolumes)
start = timer()
for i in range(0,5000):
    Numbax=Numbaenergy3D(P,bonds,orientedboundarytris,bidxTotidx,tetras,r0_ij,khook,kbend,theta0,B,MatNon,TargetVolumes)
end = timer()
print(end-start)
## Energy Minimization
start = timer()
opt.minimize(Numbaenergy3D, P.ravel()
                            ,options={'gtol':1e-01,'disp': True}  
                            ,args=(bonds
                                  ,orientedboundarytris
                                  ,bidxTotidx
                                  ,tetras
                                  ,r0_ij
                                  ,khook
                                  ,kbend
                                  ,theta0
                                  ,B
                                  ,MatNon
                                  ,TargetVolumes)
                           ).x.reshape((-1, 3))
end = timer()
print(end-start)
start = timer()
opt.minimize(energy3D, P.ravel()
                            ,options={'gtol':1e-01,'disp': True}  
                            ,args=(bonds
                                  ,orientedboundarytris
                                  ,bidxTotidx
                                  ,tetras
                                  ,r0_ij
                                  ,khook
                                  ,kbend
                                  ,theta0
                                  ,B
                                  ,MatNon
                                  ,TargetVolumes)
                           ).x.reshape((-1, 3))
end = timer()
print(end-start)
x=NumbaNeoHookean3D(r0_ij,r0_ij,khook,MatNon).sum()  
start = timer()
for i in range(0,5000):
       x=NumbaNeoHookean3D(r0_ij,r0_ij,khook,MatNon).sum()  
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>Testing the ellipse distance function<jupyter_code>import scipy.optimize as opt
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from numba import jit
def newton(f,Df,x0,epsilon,max_iter):
    '''Approximate solution of f(x)=0 by Newton's method.

    Parameters
    ----------
    f : function
        Function for which we are searching for a solution f(x)=0.
    Df : function
        Derivative of f(x).
    x0 : number
        Initial guess for a solution f(x)=0.
    epsilon : number
        Stopping criteria is abs(f(x)) < epsilon.
    max_iter : integer
        Maximum number of iterations of Newton's method.

    Returns
    -------
    xn : number
        Implement Newton's method: compute the linear approximation
        of f(x) at xn and find x intercept by the formula
            x = xn - f(xn)/Df(xn)
        Continue until abs(f(xn)) < epsilon and return xn.
        If Df(xn) == 0, return None. If the number of iterations
        exceeds max_iter, then return None.

    Examples
    --------
    >>> f = lambda x: x**2 - x - 1
    >>> Df = lambda x: 2*x - 1
    >>> newton(f,Df,1,1e-8,10)
    Found solution after 5 iterations.
    1.618033988749989
    '''
    xn = x0
    for n in range(0,max_iter):
        fxn = f(xn)
        if abs(fxn) < epsilon:
            print('Found solution after',n,'iterations.')
            return xn
        Dfxn = Df(xn)
        if Dfxn == 0:
            print('Zero derivative. No solution found.')
            return None
        xn = xn - fxn/Dfxn
    print('Exceeded maximum iterations. No solution found.')
    return None
alpha=1
beta=1
R=1
r0=np.array([0,0.5])
f=lambda theta: (alpha**2-beta**2)*R*np.sin(theta)*np.cos(theta)- alpha*r0[0]*np.sin(theta)+beta*r0[1]*np.cos(theta)
Df = lambda theta:  (alpha**2-beta**2)*R*(np.cos(theta)**2-np.sin(theta)**2)- alpha*r0[0]*np.cos(theta)-beta*r0[1]*np.sin(theta)
x0=3
epsilon=0.01
max_iter=5
newton(f,Df,x0,epsilon,max_iter)

@jit(nopython=True)
def f(theta,r0,R,alpha,beta):
    return (alpha**2-beta**2)*R*np.sin(theta)*np.cos(theta)- alpha*r0[0]*np.sin(theta)+beta*r0[1]*np.cos(theta)
@jit(nopython=True)
def Df(theta,r0,R,alpha,beta):
    return (alpha**2-beta**2)*R*(np.cos(theta)**2-np.sin(theta)**2)- alpha*r0[0]*np.cos(theta)-beta*r0[1]*np.sin(theta)

@jit(nopython=True)
def DistanceToEllipse(r0,R,alpha,beta):
    
    # Initial guess
    theta0=np.arctan2((alpha*r0[1]),(beta*r0[0]))

    # run newtons method
    max_iter=5
    theta = theta0
    for n in range(0,max_iter):
        fxn = f(theta,r0,R,alpha,beta)
        Dfxn = Df(theta,r0,R,alpha,beta)
        theta = theta - fxn/Dfxn
    
    thetafinal=theta 
    
    xellipse=R*alpha*np.cos(thetafinal)
    yellipse=R*beta*np.sin(thetafinal)
    
    deltax= r0[0]-xellipse
    deltay= r0[1]-yellipse
    
    return (thetafinal,xellipse,yellipse,np.sqrt(deltax**2+deltay**2))
start = timer()
for i in range(0,5000):
       x=DistanceToEllipse(r0,R,alpha,beta)
end = timer()
print(end-start)<jupyter_output><empty_output><jupyter_text>try out some different r0 values below. It seems to work okay!<jupyter_code>R=1
alpha=1.3
beta=1.5
theta = np.linspace(0.0, 2.0 * np.pi, 100)
x = R*alpha*np.cos(theta)
y = R*beta*np.sin(theta)
plt.plot(x,y)
plt.plot(0,0,'go')

r0=np.array([0.6,0.4])
plt.plot(r0[0],r0[1],'ro')

(thetafinal, Ellipsex,Ellipsey, distance)=DistanceToEllipse(r0,R,alpha,beta)

plt.plot(Ellipsex,Ellipsey,'bo')

# draw a ray normal to the ellipse point:
vx=-beta*np.cos(thetafinal)
vy=-alpha*np.sin(thetafinal)

plt.plot([Ellipsex,Ellipsex+vx],[Ellipsey,Ellipsey+vy])

plt.axes().set_aspect('equal')

print(distance)<jupyter_output><empty_output><jupyter_text>Testing the ellipsoid fitting functions<jupyter_code>DataFolder
# Make the Mesh
with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    #ellipsoid = geom.add_ball([0.0, 0.0, 0.0], 1)
    ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [0.95, 0.95, 1.0556])
    InputMesh = geom.generate_mesh()
    
InputMesh.write(DataFolder+"/"+"InitialMesh.vtk",binary=True) 

interiorbonds,edgebonds,boundarytris, bidxTotidx, tetras= MakeMeshData3D(InputMesh)
bonds=np.concatenate((interiorbonds,edgebonds))
orientedboundarytris=OrientTriangles(InputMesh.points,boundarytris,np.array([0,0,0]))

# Get the points on the boundary:
BoundaryPoints= np.unique(edgebonds.ravel())
P=InputMesh.points[BoundaryPoints]
xx=P[:,0]
yy=P[:,1]
zz=P[:,2]

center,axes,ec,inve,vec = ls_ellipsoid(xx,yy,zz)
print(center)
print(axes)

print(ec)
print(inve)
print(vec)<jupyter_output><empty_output><jupyter_text>This seems to work pretty well!  Testing the Analytic Prediction Functions: To test these functions are working fine, I will compare a few values to the mathematica code, with some randomly chosen parameter values:<jupyter_code>FTot<jupyter_output><empty_output><jupyter_text>Expectation:alpha0=0.5kappa0=0.3gamma0=0.2FTot(1,alpha0,gamma0=0.2,kappa0)= 16.3363FTot(1.2,alpha0,gamma0=0.2,kappa0)= 16.8611FTot(0.2,alpha0,gamma0=0.2,kappa0)= 205.077alpha0=0.5kappa0=0.3gamma0=-0.2FTot(1,alpha0,gamma0=0.2,kappa0)= 11.3097FTot(1.2,alpha0,gamma0=0.2,kappa0)= 11.7705FTot(0.2,alpha0,gamma0=0.2,kappa0)= 192.198<jupyter_code>lam=np.array([1.0001,1.2,0.2])
alpha0=0.5
kappa0=0.3
gamma0=0.2
FTot(lam,alpha0,gamma0,kappa0)
alpha0=0.5
kappa0=0.3
gamma0=-0.2
FTot(lam,alpha0,gamma0,kappa0)<jupyter_output><empty_output><jupyter_text>Conclusion: This function seems to be working okay<jupyter_code>FindGlobalMinimum<jupyter_output><empty_output><jupyter_text>Expectation:alpha0=0.6kappa0=0.4gamma0=-1FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)-> lam =1 , F =3.76991alpha0=0.6kappa0=0.4gamma0=-4.0FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)-> lam =1.55262 , F =-34.1081alpha0=0.1kappa0=0.4gamma0=-4.0FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)-> lam =0.691294 , F =-34.0731<jupyter_code>alpha0=0.6
kappa0=0.4
gamma0=-1
FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)
alpha0=0.6
kappa0=0.4
gamma0=-4.0
FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)
alpha0=0.1
kappa0=0.4
gamma0=-4.0
FindGlobalMinimum(alpha0,gamma0,kappa0,0.1,5)<jupyter_output>/home/jackbinysh/Code/ActiveElastocapillarity/Python/EnergyMinimization/AnalysisFunctions.py:175: RuntimeWarning: invalid value encountered in power
  Area = (2*np.pi/lam)*( 1+(lam**(3/2)/e)*np.arcsin(e) )
/home/jackbinysh/Code/ActiveElastocapillarity/Python/EnergyMinimization/AnalysisFunctions.py:175: RuntimeWarning: invalid value encountered in arcsin
  Area = (2*np.pi/lam)*( 1+(lam**(3/2)/e)*np.arcsin(e) )
/home/jackbinysh/Code/ActiveElastocapillarity/Python/EnergyMinimization/AnalysisFunctions.py:177: RuntimeWarning: invalid value encountered in arctanh
  Fbend=(2/3)*np.pi*kappa*(7+(2/lam**3)+3*lam**3*np.arctanh(np.lib.scimath.sqrt(1-lam**3))/np.lib.scimath.sqrt(1-lam**3))
{'fraction_non_alphanumeric': 0.1447574730666127, 'fraction_numerical': 0.04612816751808204, 'mean_word_length': 4.393344244408074, 'pattern_counts': {'":': 0, '<': 95, '<?xml version=': 0, '>': 104, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '665967', 'license': 'MIT', 'n_tokens_mistral': 8313, 'n_tokens_neox': 7794, 'n_words': 1281, 'path': 'Python/EnergyMinimization/UnitTests.ipynb', 'repo_name': 'SouslovLab/ActiveElastocapillarity', 'size': '19771'}