"""Original grouped Bernoulli mixture, two actual E/M updates and exact likelihood oracle."""
from fractions import Fraction as F
from itertools import product
import math,json

def validate(groups,pi,ps):
 if type(groups) is not list or not 1<=len(groups)<=6 or any(type(g) is not list or not 1<=len(g)<=4 or any(type(x) is not int or x not in [0,1] for x in g) for g in groups) or len({len(g) for g in groups})!=1:raise ValueError('1..6 equal-length binary groups,1..4 flips each')
 if any(type(a) is not list or len(a)!=2 or any(type(x) not in [int,F] for x in a) for a in [pi,ps]) or sum(pi)!=1 or any(x<=0 for x in pi) or any(not 0<=x<=1 for x in ps):raise ValueError('two positive normalized mixture masses and Bernoulli parameters0..1')

def contribution(group,pi,ps,k):
 h=sum(group);m=len(group);return F(pi[k])*F(ps[k])**h*(1-F(ps[k]))**(m-h)

def likelihood(groups,pi,ps):
 validate(groups,pi,ps);v=F(1)
 for g in groups:
  mass=sum(contribution(g,pi,ps,k) for k in range(2))
  if mass==0:raise ValueError('zero observation likelihood')
  v*=mass
 return v

def latent_oracle(groups,pi,ps):
 total=F(0)
 for assignment in product(range(2),repeat=len(groups)):
  mass=F(1)
  for g,k in zip(groups,assignment):mass*=contribution(g,pi,ps,k)
  total+=mass
 return total

def step(groups,pi,ps):
 before=likelihood(groups,pi,ps);weights=[]
 for g in groups:
  vals=[contribution(g,pi,ps,k) for k in range(2)];r=[v/sum(vals) for v in vals];assert sum(r)==1;weights.append(r)
 n=len(groups);m=len(groups[0]);nk=[sum(row[k] for row in weights) for k in range(2)]
 if any(z==0 for z in nk):raise ValueError('zero component responsibility under local contract')
 heads=[sum(weights[i][k]*sum(g) for i,g in enumerate(groups)) for k in range(2)];newpi=[z/n for z in nk];newps=[heads[k]/(m*nk[k]) for k in range(2)];after=likelihood(groups,newpi,newps)
 assert after>=before and after==latent_oracle(groups,newpi,newps)
 # Weighted Bernoulli score equation zero and strictly negative second derivative.
 scores=[]
 for k,p in enumerate(newps):
  tails=m*nk[k]-heads[k]
  if 0<p<1:
   derivative=heads[k]/p-tails/(1-p);second=-heads[k]/p**2-tails/(1-p)**2;assert derivative==0 and second<0;scores.append({'component':k,'score':str(derivative),'second_derivative':str(second)})
 encode=lambda xs:[str(x) for x in xs]
 row={'old_pi':encode(pi),'old_p':encode(ps),'responsibility_rows':[encode(r) for r in weights],'Nk':encode(nk),'weighted_heads':encode(heads),'correct_denominators_m_times_Nk':encode([m*x for x in nk]),'new_pi':encode(newpi),'new_p':encode(newps),'wrong_divide_heads_by_number_of_groups':encode([x/n for x in heads]),'old_exact_likelihood':str(before),'new_exact_likelihood':str(after),'old_log_likelihood':math.log(before.numerator)-math.log(before.denominator),'new_log_likelihood':math.log(after.numerator)-math.log(after.denominator),'independent_all_latent_assignments':2**n,'weighted_Bernoulli_score_oracle':scores}
 return newpi,newps,row

def demo():
 groups=[[0,0,0],[0,0,1],[1,1,0],[1,1,1]];pi=[F(1,2)]*2;ps=[F(1,4),F(3,4)];initial=likelihood(groups,pi,ps);assert initial==latent_oracle(groups,pi,ps);rows=[]
 for i in range(2):pi,ps,row=step(groups,pi,ps);rows.append(row)
 assert rows[0]['new_p']==['19/84','65/84'] and rows[1]['new_p']==['60955/281484','220529/281484']
 spi,sps,sym=step(groups,[F(1,2)]*2,[F(1,2)]*2);assert sps==[F(1,2)]*2 and sym['old_exact_likelihood']==sym['new_exact_likelihood']
 swap_pi,swap_ps,swap=step(groups,[F(1,2)]*2,[F(3,4),F(1,4)]);assert swap_ps==list(reversed(list(map(F,rows[0]['new_p']))))
 rejected=[]
 for label,g,pp,qq in [('empty',[],[F(1,2)]*2,[F(1,4),F(3,4)]),('nonbinary',[[0,2]],[F(1,2)]*2,[F(1,4),F(3,4)]),('unequal_lengths',[[0],[1,0]],[F(1,2)]*2,[F(1,4),F(3,4)]),('invalid_pi',groups,[F(1,4)]*2,[F(1,4),F(3,4)]),('invalid_p',groups,[F(1,2)]*2,[-1,2]),('zero_likelihood',[[1]],[F(1,2)]*2,[0,0])]:
  try:step(g,pp,qq)
  except ValueError:rejected.append(label)
  else:raise AssertionError('expected rejection')
 tiny=F(1,10**400);tpi,tps,trow=step([[1],[1]],[F(1,2)]*2,[tiny,tiny]);assert F(trow['old_exact_likelihood'])>0 and float(F(trow['old_exact_likelihood']))==0 and math.isfinite(trow['old_log_likelihood']) and abs(trow['old_log_likelihood']+800*math.log(10))<1e-9 and trow['new_log_likelihood']==0
 return {'positive_likelihood_float_underflow_log_regression':trow,'log_display_computed_from_integer_numerator_minus_denominator':True,'groups':groups,'latent_contract':'one hidden component per complete three-flip group, not per individual flip','initial_exact_likelihood':str(initial),'two_actual_EM_steps':rows,'symmetric_initialization_stays_symmetric':sym,'label_swap_step':swap,'invalid_cases_rejected':rejected,'scope':'observed ordered binary groups, no binomial multiplicity factor; Fraction responsibilities/parameters/likelihood and floating natural log display; two likelihood improvements not global optimum/true latent coin discovery/real model training'}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=True,indent=2))