1import json, math, random
 2import numpy as np
 3import torch
 4from scipy.linalg import eigh
 5
 6SEED=2950
 7np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED)
 8
 9def beam(nel=8, L=1.0, EI=1.0, rhoA=1.0):
10    h=L/nel; K=np.zeros((2*(nel+1),2*(nel+1))); M=K.copy()
11    ke=EI/h**3*np.array([[12,6*h,-12,6*h],[6*h,4*h*h,-6*h,2*h*h],[-12,-6*h,12,-6*h],[6*h,2*h*h,-6*h,4*h*h]],float)
12    me=rhoA*h/420*np.array([[156,22*h,54,-13*h],[22*h,4*h*h,13*h,-3*h*h],[54,13*h,156,-22*h],[-13*h,-3*h*h,-22*h,4*h*h]],float)
13    for e in range(nel):
14        ix=[2*e,2*e+1,2*e+2,2*e+3]; K[np.ix_(ix,ix)]+=ke; M[np.ix_(ix,ix)]+=me
15    # clamped left, free right
16    return K[2:,2:], M[2:,2:]
17
18def inertia(K,M,w,tol=1e-12):
19    # For positive-definite M, inertia(K-w^2 M) equals count of frequencies below w.
20    A=K-(w*w)*M; ev=np.linalg.eigvalsh(A); scale=max(1.,np.max(np.abs(A)))
21    return int(np.sum(ev < -tol*scale))
22
23def bracket(K,M,target,hi=None,tol=1e-3):
24    vals=eigh(K,M,eigvals_only=True); freqs=np.sqrt(np.maximum(vals,0))
25    if hi is None: hi=float(freqs[-1]*1.05)
26    lo=0.; upper=hi
27    # Find the transition count target-1 -> target; bisection around exact transition.
28    grid=np.linspace(0,upper,400)
29    for a,b in zip(grid[:-1],grid[1:]):
30        if inertia(K,M,a)==target-1 and inertia(K,M,b)>=target:
31            lo,upper=float(a),float(b); break
32    # narrow an interval whose endpoints have counts target-1 and target
33    while upper-lo>tol:
34        mid=(lo+upper)/2
35        if inertia(K,M,mid)<target: lo=mid
36        else: upper=mid
37    return lo,upper,freqs
38
39def train(K,M,target,interval=None,steps=2000,seed=0):
40    torch.manual_seed(seed); np.random.seed(seed)
41    dev='cuda' if torch.cuda.is_available() else 'cpu'
42    try:
43        kt=torch.tensor(K,dtype=torch.float64,device=dev); mt=torch.tensor(M,dtype=torch.float64,device=dev)
44        n=K.shape[0]; q=torch.randn(n,dtype=torch.float64,device=dev,requires_grad=True)
45        # omega starts randomly in a broad range; bracketed version starts inside its certified interval
46        exact=np.sqrt(eigh(K,M,eigvals_only=True));
47        if interval is None: w=torch.tensor(float(np.random.uniform(.2,exact[-1]*1.15)),dtype=torch.float64,device=dev,requires_grad=True)
48        else: w=torch.tensor(float((interval[0]+interval[1])/2),dtype=torch.float64,device=dev,requires_grad=True)
49        opt=torch.optim.Adam([q,w],lr=.025)
50        for t in range(steps):
51            nq=q/(torch.linalg.vector_norm(q)+1e-12)
52            r=kt@nq-(w*w)*(mt@nq)
53            loss=torch.sum(r*r)
54            # This is the stated box penalty; modest normalization avoids trivial q=0.
55            if interval is not None:
56                lo,up=interval; loss=loss+2000*(torch.relu(torch.tensor(lo,device=dev)-w)**2+torch.relu(w-torch.tensor(up,device=dev))**2)
57                if t < steps//5: # freeze frequency during initial eigenvector settling
58                    w.grad=None
59            loss.backward(); opt.step(); opt.zero_grad()
60            if interval is not None and t < steps//5: w.data.fill_((interval[0]+interval[1])/2)
61        qn=(q.detach()/torch.linalg.vector_norm(q.detach())).cpu().numpy(); ww=float(w.detach().cpu())
62        ev,vec=eigh(K,M); freqs=np.sqrt(ev); # mass-normalized reference vectors
63        mac=np.array([(abs(qn@vec[:,i])**2/((qn@qn)*(vec[:,i]@vec[:,i]))) for i in range(len(freqs))])
64        return ww, int(np.argmax(mac)), float(mac.max()), float(freqs[np.argmin(abs(freqs-ww))])
65    except Exception:
66        # CPU retry is required for shared/fragile CUDA environments.
67        if dev=='cuda':
68            torch.cuda.empty_cache(); torch.cuda.is_available=lambda: False
69            return train(K,M,target,interval,steps,seed)
70        raise
71
72def main():
73    K,M=beam(8); exact=np.sqrt(eigh(K,M,eigvals_only=True)); target=3
74    ws=np.linspace(0,exact[-1]*1.1,250)
75    counts=np.array([inertia(K,M,float(w)) for w in ws])
76    monotone=bool(np.all(np.diff(counts)>=0)); changes=int(np.sum(np.diff(counts)>0))
77    lo,up,freqs=bracket(K,M,target)
78    contained=int(np.sum((freqs>lo)&(freqs<up)))
79    baseline=[train(K,M,target,None,2000,s) for s in range(8)]
80    idea=[train(K,M,target,(lo,up),2000,s) for s in range(8)]
81    # lower-mode collapse means identified FEM mode index below requested target.
82    b_lower=sum(x[1]<target-1 for x in baseline); i_lower=sum(x[1]<target-1 for x in idea)
83    out={'exact_frequencies':exact.tolist(),'count_monotone':monotone,'count_grid_changes':changes,'bracket':[lo,up],'bracket_width':up-lo,'frequencies_inside_bracket':contained,'target_mode_zero_based':target-1,
84         'baseline':baseline,'idea':idea,'baseline_lower_mode_rate':b_lower/len(baseline),'idea_lower_mode_rate':i_lower/len(idea),
85         'baseline_target_rate':sum(x[1]==target-1 for x in baseline)/len(baseline),'idea_target_rate':sum(x[1]==target-1 for x in idea)/len(idea)}
86    print(json.dumps(out,indent=2))
87if __name__=='__main__': main()