"""Scientific tests of source conditioning, generator, budgets and path accounts."""
import unittest
from fractions import Fraction as F
from collections import Counter
from itertools import combinations
import math
import numpy as np
import mpmath as mp
from polymer import *
from source import *
from certificates import *


class ScientificChecks(unittest.TestCase):
    def test_01_catalogue_census_and_literal_generator(self):
        for n in [4,5,6]:
            cat=Catalogue(n);self.assertEqual(len(cat.words),2**(n+1)-2);self.assertEqual(len(cat.splits),(n-2)*2**(n+1)+4)
            labels=cat.productive_labels();self.assertEqual(len(labels),224);self.assertEqual(sum(z in FOOD for z,r in labels),192)
        cat=Catalogue();env=Environment.matched(cat);reactor=PolymerReactor(cat,env)
        self.assertEqual(len(reactor.channels),174)
        for ch in reactor.channels:
            if ch.kind=='internal':self.assertEqual(sum(len(z)*d for z,d in ch.change.items()),0)
        repeated=Channel(('00','00','00'),('0000','00'),4.)
        self.assertAlmostEqual(repeated.propensity({'00':3},10),.24,14)
        self.assertEqual(repeated.propensity({'00':2},10),0)
        reverse=Channel(('0011','0011'),('00','11','0011'),4.)
        self.assertEqual(reverse.propensity({'0011':1},10),0)
        counts=reactor.initial(100);counts[cat.index['0011']]=7
        vector=reactor.count_rates(counts,100);mapping=dict(zip(cat.words,counts))
        for value,ch in zip(vector,reactor.channels):self.assertAlmostEqual(value,ch.propensity(mapping,100),12)
        self.assertTrue(env.witness_class(cat));self.assertFalse(Environment.matched(cat,'food').witness_class(cat))

    def test_02_capped_law_conditional_likelihood_and_shared_marks(self):
        source=CappedZipf('1.5');R=6;p=source.degree_probabilities(R)
        with mp.workdps(55):
            self.assertLess(abs(sum(p)-1),mp.mpf('1e-50'))
            conditional,evidence=source.conditional_degrees(R,1,1)
            brute=[mp.mpf(0)]*R
            for d in range(R):
                subsets=list(combinations(range(R),d))
                brute[d]=p[d]*sum(0 in s and 1 not in s for s in subsets)/len(subsets)
            total=sum(brute)
            for a,b in zip(conditional,brute):self.assertLess(abs(a-b/total),mp.mpf('1e-50'))
        cat=Catalogue();env=source.sample_environment(cat,59123,True);env.validate(cat);self.assertTrue(env.witness_class(cat))
        self.assertFalse(any(z in FOOD for z,r,H in env.incidences))
        reactor=PolymerReactor(cat,env)
        # Every basal/catalytic reversible pair shares exactly one coefficient.
        internal=[ch for ch in reactor.channels if ch.kind=='internal']
        for a,b in zip(internal[::2],internal[1::2]):self.assertEqual(a.coefficient,b.coefficient);self.assertEqual(a.inputs,b.outputs)
        self.assertEqual(env.deleted().S,env.S);self.assertEqual(env.deleted().B,env.B)
        row=source.sample_row(68,np.random.default_rng(5),required=(0,1),forbidden=(2,3));self.assertTrue({0,1}.issubset(row));self.assertFalse({2,3}&set(row))
        with self.assertRaises(ValueError):source.sample_row(68,np.random.default_rng(5),(1,),(1,))

    def test_03_moments_pair_dependence_and_finite_converse(self):
        with mp.workdps(60):
            for a in ['1.5','2','2.5','4']:
                source=CappedZipf(a);m=source.moments(4);p=source.degree_probabilities(68)
                ED=sum(d*v for d,v in enumerate(p));ED2=sum(d*(d-1)*v for d,v in enumerate(p))
                self.assertLess(abs(m['mean_degree']-ED),mp.mpf('1e-45'));self.assertLess(abs(m['second_factorial_moment']-ED2),mp.mpf('1e-43'))
                self.assertGreater(m['same_row_pair'],m['incidence']**2)
            m=CappedZipf('1.5').moments(4)
            self.assertGreater(float(m['conditional_cap_atom']),.5);self.assertLess(float(m['cap_atom']),.1)
            heavy=CappedZipf('1.5').moments(64);light=CappedZipf('2.5').moments(64)
            self.assertAlmostEqual(float(heavy['pair_ratio']),2/3,6);self.assertLess(float(light['pair_ratio']),1e-8)
            residual=CappedZipf('1.5').small_cap_residual(4)
            self.assertAlmostEqual(float(residual['log10_complement']),-11.446,3)
            n=15;R=(n-2)*2**(n+1)+4;a=mp.mpf('1.5')
            direct=(mp.zeta(a-1)-mp.zeta(a-1,R)+R*mp.zeta(a,R))/mp.zeta(a)-1
            self.assertLess(abs(CappedZipf('1.5').moments(n)['mean_degree']/direct-1),mp.mpf('1e-45'))
        self.assertEqual(screening_trials('.001'),2995)
        # Fresh environments and repeats in one environment differ even at the
        # same average success: q in {0,1}, equal mass, N=2 gives .5 versus .75.
        self.assertNotEqual(1-(.5*1**2+.5*0**2),1-(1-.5)**2)

    def test_04_startup_necessary_scale_and_corrected_scalar_floor(self):
        c=deterministic_noise_checks();self.assertEqual(c['mass_upper_from_prefix_noise'],F(21,2))
        self.assertGreater(c['product_at_one_lower'],c['stock_floor']);self.assertGreater(c['product_forcing_lower'],F(1,10**15))
        self.assertEqual(c['interval99_export_lower'],F(319062203,2*10**9))
        self.assertLess(c['interval78_export_lower'],F(1,10));self.assertGreater(c['interval79_export_lower'],F(1,10))
        for y in map(F,[0,1,2,F(11,3),F(11,2)]):self.assertEqual(F(121,3)-(22*y-3*y*y),3*(y-F(11,3))**2)
        self.assertEqual(necessary_startup_scale()['necessary_integer_scale'],14272222)
        for V in [100,200]:
            b=first_birth_ceiling(V);value=-math.expm1(-(484/3)*2e-9*V)
            self.assertLessEqual(float(b['lower_evaluation'])-1e-15,value);self.assertGreaterEqual(float(b['upper_evaluation'])+1e-15,value)
        # Demonstrate why the compensator is not fixed even when mean mass is 10.
        # L=10+eta(1-exp(-t)) has prefix noise eta*t, not constant compensator.
        eta=.01;t=.5;L=10+eta*(1-math.exp(-t));compensator=10-eta*(t+math.exp(-t)-1)
        self.assertNotEqual(compensator,10);self.assertAlmostEqual(L-compensator,eta*t,13)

    def test_05_six_term_enclosure_and_scale_hypotheses(self):
        c=TwoWindowCertificate();result=c.evaluate();self.assertEqual(result['status'],'valid conditional bound');self.assertLess(result['failure_upper'],F(1,100))
        scale=c.sufficient_scale();self.assertEqual(scale['dominant_term'],'selected_product');self.assertTrue(5e49<float(scale['count_scale'])<5.001e49)
        self.assertIsNone(TwoWindowCertificate(V=100).evaluate()['failure_upper'])
        self.assertFalse(c.uniform()['envelope_hypothesis'])
        big=TwoWindowCertificate(V=10**60*25);self.assertTrue(big.uniform()['envelope_hypothesis']);self.assertGreater(big.uniform()['exponent'],48)
        # The appendix's equality to 48 is not used: 48 is only a lower bound.
        self.assertNotEqual(big.uniform()['exponent'],48)
        with self.assertRaises(ValueError):TwoWindowCertificate(tolerance_fraction=F(1))

    def test_06_matched_ode_and_signed_material_ledger(self):
        cat=Catalogue();times=np.array([0.,1.,100.,199.]);X=len(cat.words);selected=cat.index['0011']
        rows={}
        for kind in ['self','food','deleted']:
            reactor=PolymerReactor(cat,Environment.matched(cat,kind));y=reactor.deterministic(times);rows[kind]=y
            self.assertLess(max(abs(y[:,:X]@reactor.length-10)),1e-9)
            self.assertLess(max(abs(y[:,:X]@reactor.nonfood+y[:,X]-y[:,X+3])),1e-9)
            self.assertAlmostEqual(y[-1,X+1],1194,7);self.assertAlmostEqual(y[-1,X+2],1990,7)
        self.assertAlmostEqual(rows['self'][-1,selected],.2753,4)
        self.assertAlmostEqual(rows['food'][-1,selected],.3441,4)
        self.assertTrue(all(v>100 for v in np.diff(rows['self'][1:,X])))
        self.assertGreater(rows['deleted'][-1,X],0);self.assertLess(rows['deleted'][-1,X],.0001)

    def test_07_discrete_path_budget_and_whole_mission(self):
        cat=Catalogue();reactor=PolymerReactor(cat,Environment.matched(cat))
        short=reactor.simulate(20,1,5901,10000);self.assertTrue(short['completed']);self.assertEqual(short['ledger_residual'],0);self.assertIsNone(short['mission'])
        incomplete=reactor.simulate(20,199,5901,1);self.assertFalse(incomplete['completed']);self.assertIsNone(incomplete['mission'])
        food=PolymerReactor(cat,Environment.matched(cat,'food')).simulate(100,.1,5902,10000)
        self.assertTrue(food['completed']);self.assertIsNotNone(food['first_nonfood']);self.assertEqual(food['first_nonfood'],food['first_selected'])
        self.assertEqual(food['final_nonfood']+food['export'],food['signed_synthesis'])
        self.assertEqual(F(1,5)/(10+2*1195),F(1,12000))


if __name__=='__main__':unittest.main()