"""Original small Euclidean NMF update exercise, existing NumPy only."""
import numpy as np
import json,math
def matrix(a):
 if type(a) is not list or not 1<=len(a)<=4 or any(type(r) is not list or not 1<=len(r)<=4 for r in a) or len({len(r) for r in a})!=1:raise ValueError('bounded rectangular list')
 if any(type(x) not in [int,float] or abs(x)>1000 or not math.isfinite(x) or x<0 for r in a for x in r):raise ValueError('finite nonnegative bounded reals')
 return np.array(a,dtype=np.float64)
def multiply_loop(a,b):return [[sum(float(a[i,k])*float(b[k,j]) for k in range(a.shape[1])) for j in range(b.shape[1])] for i in range(a.shape[0])]
def step(v,w,h):
 v=matrix(v);w=matrix(w);h=matrix(h)
 if v.shape!=(w.shape[0],h.shape[1]) or w.shape[1]!=h.shape[0]:raise ValueError('shape mismatch')
 hn=w.T@v;hd=w.T@w@h
 if np.any(hd<=0):raise ValueError('zero denominator H, no epsilon added')
 hh=h*(hn/hd);wn=v@hh.T;wd=w@hh@hh.T
 if np.any(wd<=0):raise ValueError('zero denominator W, no epsilon added')
 ww=w*(wn/wd)
 assert np.allclose(multiply_loop(w,h),w@h,atol=1e-12,rtol=0)
 assert np.allclose(multiply_loop(ww,hh),ww@hh,atol=1e-12,rtol=0)
 return ww,hh,{'H_numerator':hn.tolist(),'H_denominator':hd.tolist(),'H_after':hh.tolist(),'W_numerator_using_updated_H':wn.tolist(),'W_denominator_using_updated_H':wd.tolist()}
def loss(v,w,h):
 a=np.array(v);ww=np.array(w);hh=np.array(h);computed=float(np.sum((a-ww@hh)**2));ref=multiply_loop(ww,hh)
 independent=sum((float(a[i,j])-ref[i][j])**2 for i in range(a.shape[0]) for j in range(a.shape[1]))
 assert abs(computed-independent)<=1e-12
 return computed
def demo():
 v=[[3,1,0],[1,3,1],[0,1,3]];w=[[1,.5],[.5,1],[1,1]];h=[[1,.5,.8],[.5,1,.2]];initial={'V':v,'W':w,'H':h};history=[];old=loss(v,w,h)
 for i in range(5):
  ww,hh,details=step(v,w,h);half=loss(v,w,hh);new=loss(v,ww,hh);assert half<=old+1e-12 and new<=half+1e-12
  history.append({'iteration':i+1,'loss_before':old,'loss_after_H':half,'loss_after_W':new,'W_after':ww.tolist(),**details});w=ww.tolist();h=hh.tolist();old=new
 d=np.array([2,.5]);scaled_w=np.array(w)*d;scaled_h=np.array(h)/d[:,None];assert np.allclose(scaled_w@scaled_h,np.array(w)@np.array(h),atol=1e-12,rtol=0)
 zw=[r[:] for r in initial['W']];zw[0][0]=0;zww,zhh,_=step(v,zw,initial['H']);assert zww[0,0]==0
 rejected=[]
 bad=[('negative',[[-1]],[[1]],[[1]]),('nan',[[float('nan')]],[[1]],[[1]]),('shape',[[1,2]],[[1]],[[1]]),('zero_denominator',[[1]],[[0]],[[1]])]
 for label,a,b,c in bad:
  try:step(a,b,c)
  except ValueError:rejected.append(label)
  else:raise AssertionError('expected rejection')
 return {'numpy_version':np.__version__,'initial':initial,'five_actual_H_then_W_updates':history,'loss_definition':'sum of squared entries V-WH, no normalization','nonnegative_factors':True,'loop_matmul_oracles_equal_at_1e_minus12_absolute_tolerance':True,'scale_D':[2,.5],'scaled_W':scaled_w.tolist(),'inverse_scaled_H':scaled_h.tolist(),'scaled_reconstruction_equal':True,'zero_lock_W00_after_one_update':float(zww[0,0]),'rejected':rejected,'limits':'no epsilon, no global optimum, no topic quality or trained product/model claim'}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=True,indent=2))