verify_eqangles.py

import math
import numpy as np
import itertools
from typing import List, Tuple, Optional
from relations import predicate

# Types come from your codebase
# - GeometricProblem has .derived_facts (list of predicate objects)
# - Points have .name, .x, .y
# - Eqangle predicate contains 6 points representing two angles (A,O,B,C,O',D)


def _get_points_from_pred(pred) -> Optional[List]:
    """Best-effort to extract point list from a predicate.
    Returns list of points (possibly 6 for Eqangle) or None.
    """
    # Preferred unified API if available
    if hasattr(pred, 'get_points'):
        try:
            pts = pred.get_points()
            if isinstance(pts, list) and pts and hasattr(pts[0], 'x'):
                return pts
        except Exception:
            pass
    # Fallback common attr
    if hasattr(pred, 'points'):
        pts = getattr(pred, 'points')
        if isinstance(pts, (list, tuple)) and pts and hasattr(pts[0], 'x'):
            return list(pts)
    return None


def _angle_from_points(A, O, B) -> float:
    """Compute directed angle AOB in radians in [0, pi] using the diagram.
    Mirrors the logic used in ar.py._get_angle_from_diagram.
    Returns -1 if degenerate (zero vector).
    """
    OA = np.array([A.x - O.x, A.y - O.y], dtype=float)
    OB = np.array([B.x - O.x, B.y - O.y], dtype=float)
    na = np.linalg.norm(OA)
    nb = np.linalg.norm(OB)
    if na == 0.0 or nb == 0.0:
        return -1.0
    det = OA[0] * OB[1] - OA[1] * OB[0]
    dot = float(np.dot(OA, OB))
    cos_val = max(-1.0, min(1.0, dot / (na * nb)))
    theta = math.acos(cos_val)
    return theta if det > 0 else math.pi - theta
def _equal_mod_pi(a: float, b: float, tol: float = 1e-4) -> bool:
    if a < 0 or b < 0:
        return False
    diff = abs(a - b)
    diff = min(diff, math.pi - diff)
    return diff <= tol

def _format_angle_name(A, O, B) -> str:
    return f"{A.name}{O.name}{B.name}"

def _segment_length(P, Q) -> float:
    return float(math.hypot(P.x - Q.x, P.y - Q.y))

def _seg_endpoints(seg):
    try:
        if isinstance(seg, (list, tuple)) and len(seg) == 2:
            return seg[0], seg[1]
        pts = list(seg)
        if len(pts) != 2:
            return None, None
        pts = sorted(pts, key=lambda p: p.name)
        return pts[0], pts[1]
    except Exception:
        return None, None

def _collinear(A, B, C, tol: float = 1e-8) -> bool:
    ABx, ABy = (B.x - A.x), (B.y - A.y)
    ACx, ACy = (C.x - A.x), (C.y - A.y)
    cross = ABx * ACy - ABy * ACx
    return abs(cross) <= tol

def _parallel(A, B, C, D, tol: float = 1e-8) -> bool:
    v1x, v1y = (B.x - A.x), (B.y - A.y)
    v2x, v2y = (D.x - C.x), (D.y - C.y)
    cross = v1x * v2y - v1y * v2x
    return abs(cross) <= tol

def _sameclock(a1, o1, b1, a2, o2, b2) -> bool:
    def orient(a, o, b):
        return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
    s1 = orient(a1, o1, b1)
    s2 = orient(a2, o2, b2)
    if s1 == 0 or s2 == 0:
        return False
    return (s1 > 0 and s2 > 0) or (s1 < 0 and s2 < 0)

def _fit_circle_through_three(P1, P2, P3, tol: float = 1e-12):
    """Compute circle center and radius through three non-collinear points."""
    x1, y1 = P1.x, P1.y
    x2, y2 = P2.x, P2.y
    x3, y3 = P3.x, P3.y
    a = x1 * (y2 - y3) - y1 * (x2 - x3) + x2 * y3 - x3 * y2
    if abs(a) < tol:
        return None, None
    b = (x1**2 + y1**2) * (y3 - y2) + (x2**2 + y2**2) * (y1 - y3) + (x3**2 + y3**2) * (y2 - y1)
    c = (x1**2 + y1**2) * (x2 - x3) + (x2**2 + y2**2) * (x3 - x1) + (x3**2 + y3**2) * (x1 - x2)
    cx = -b / (2 * a)
    cy = -c / (2 * a)
    R = math.hypot(x1 - cx, y1 - cy)
    return (cx, cy), R

def _on_circle(center, R, P, tol: float = 1e-6) -> bool:
    """Check if point P lies on circle with center and radius."""
    cx, cy = center
    return math.isclose(math.hypot(P.x - cx, P.y - cy), R, rel_tol=tol, abs_tol=tol)


def _predicate_ok(pred, tol: float) -> bool:
    """Return True if predicate holds in the diagram, False otherwise.
    Supports: Eqangle, Eqratio, Simtri1/2, Contri1/2, Para, Perp, Cong, Midp, Sameclock, Cyclic, Circle.
    """
    pts = _get_points_from_pred(pred)
    if pts is None:
        return True
    cls_name = type(pred).__name__.lower()

    # Eqangle
    if 'eqangle' in cls_name and len(pts) == 6:
        A, O, B, C, O2, D = pts
        val1 = _angle_from_points(A, O, B)
        val2 = _angle_from_points(C, O2, D)
        return _equal_mod_pi(val1, val2, tol)

    # Eqratio
    if isinstance(pred, predicate.Eqratio):
        try:
            segs = pred.ratios
            r1, r2 = segs[0]
            r3, r4 = segs[1]
            ab = _segment_length(r1[0], r1[1])
            cd = _segment_length(r2[0], r2[1])
            ef = _segment_length(r3[0], r3[1])
            gh = _segment_length(r4[0], r4[1])
            return (cd != 0 and gh != 0 and math.isclose(ab / cd, ef / gh, rel_tol=1e-6, abs_tol=1e-6))
        except Exception:
            return False

    # Congruent triangles only (ordered)
    if isinstance(pred, predicate.Contri1):
        return _verify_congruent(pred.triangle1, pred.triangle2, ordered=True)

    # Para
    if isinstance(pred, predicate.Para) and len(pts) == 4:
        A, B, C, D = pts
        return _parallel(A, B, C, D)

    # Perp
    if 'perp' in cls_name and len(pts) == 4:
        P1, P2, P3, P4 = pts
        v1 = np.array([P2.x - P1.x, P2.y - P1.y], dtype=float)
        v2 = np.array([P4.x - P3.x, P4.y - P3.y], dtype=float)
        n1 = np.linalg.norm(v1)
        n2 = np.linalg.norm(v2)
        if n1 == 0 or n2 == 0:
            return False
        det = v1[0] * v2[1] - v1[1] * v2[0]
        dot = float(np.dot(v1, v2))
        ang = abs(math.atan2(det, dot))
        near = min(abs(ang - math.pi/2), abs((math.pi - ang) - math.pi/2))
        return near <= tol

    # Cong
    if 'cong' in cls_name and len(pts) == 4:
        A, B, C, D = pts
        return abs(_segment_length(A, B) - _segment_length(C, D)) <= (1e-6 if tol < 1e-6 else tol)

    # Midp
    if 'midp' in cls_name and len(pts) == 3:
        M, A, B = pts
        d1 = _segment_length(M, A)
        d2 = _segment_length(M, B)
        equal = abs(d1 - d2) <= (1e-6 if tol < 1e-6 else tol)
        return equal and _collinear(A, B, M)

    # Sameclock
    if isinstance(pred, predicate.Sameclock) and len(pts) == 6:
        a1, o1, b1, a2, o2, b2 = pts
        return _sameclock(a1, o1, b1, a2, o2, b2)

    # Cyclic
    if isinstance(pred, predicate.Cyclic) and len(pts) == 4:
        P1, P2, P3, P4 = pts
        center, R = _fit_circle_through_three(P1, P2, P3)
        return (center is not None and R is not None and _on_circle(center, R, P4))

    # Circle
    if isinstance(pred, predicate.Circle) and len(pts) == 4:
        C, P1, P2, P3 = pts
        R = _segment_length(C, P1)
        return (_on_circle((C.x, C.y), R, P2) and _on_circle((C.x, C.y), R, P3))

    return True


def verify_eqangles(problem, outfile_path: str = "eqangle_verification.txt", tol_deg: float = 0.01) -> int:
    """Verify all Eqangle predicates in problem.derived_facts against the diagram.
    Writes one line per eqangle:
      eqangle AOB (X.Y deg), COD (U.V deg) -> equal_mod_pi=true/false
    Returns the number of eqangle predicates checked.
    """
    if not hasattr(problem, 'derived_facts'):
        raise ValueError("Problem has no derived_facts to verify.")

    tol = math.radians(tol_deg)
    # Build indices of reasoning steps: by object equality and by string
    reason_index_obj = {}
    reason_index_str = {}
    try:
        for iter_level in getattr(problem, 'reasoning_steps', []) or []:
            for entry in iter_level or []:
                if not isinstance(entry, dict):
                    continue
                pr = entry.get('predicate')
                deps = entry.get('dependencies', {}) or {}
                if pr is not None:
                    try:
                        reason_index_obj.setdefault(pr, []).append(deps)
                    except Exception:
                        pass
                    try:
                        reason_index_str.setdefault(str(pr), []).append(deps)
                    except Exception:
                        pass
    except Exception:
        reason_index_obj, reason_index_str = {}, {}
    # Index derived facts by their string for lookup of actual objects
    pred_index_by_str = {}
    try:
        for df in getattr(problem, 'derived_facts', []) or []:
            try:
                pred_index_by_str[str(df)] = df
            except Exception:
                pass
    except Exception:
        pred_index_by_str = {}

    count = 0
    lines: List[str] = []

    for pred in getattr(problem, 'derived_facts', []):
        cls_name = type(pred).__name__.lower()
        pts = _get_points_from_pred(pred)
        if not pts:
            continue

        # Eqangle: 6 points A,O,B,C,O',D
        if 'eqangle' in cls_name and len(pts) == 6:
            A, O, B, C, O2, D = pts
            val1 = _angle_from_points(A, O, B)
            val2 = _angle_from_points(C, O2, D)
            equal = _equal_mod_pi(val1, val2, tol)
            deg1 = val1 * 180.0 / math.pi if val1 >= 0 else float('nan')
            deg2 = val2 * 180.0 / math.pi if val2 >= 0 else float('nan')
            name1 = _format_angle_name(A, O, B)
            name2 = _format_angle_name(C, O2, D)

            # Collect dependency sources for this eqangle (e.g., DD or AR)
            sources = set()
            deps_list = []
            try:
                deps_list.extend(reason_index_obj.get(pred, []) or [])
            except Exception:
                pass
            try:
                deps_list.extend(reason_index_str.get(str(pred), []) or [])
            except Exception:
                pass
            for deps in deps_list:
                src = deps.get('type') or deps.get('source') or 'unknown'
                sources.add(str(src))
            source_str = '|'.join(sorted(sources)) if sources else 'unknown'

            if not equal:
                lines.append(
                    f"eqangle {name1} ({deg1:.4f} deg), {name2} ({deg2:.4f} deg) -> equal_mod_pi={equal} | source={source_str}"
                )
            count += 1
            continue

        # Eqratio: verify AB/CD == EF/GH via lengths
        if isinstance(pred, predicate.Eqratio):
            try:
                segs = pred.ratios  # [[seg1, seg2],[seg3, seg4]], each seg is list of two Points
                r1, r2 = segs[0]
                r3, r4 = segs[1]
                A,B = _seg_endpoints(r1)
                C,D = _seg_endpoints(r2)
                E,F = _seg_endpoints(r3)
                G,H = _seg_endpoints(r4)
                if any(p is None for p in [A,B,C,D,E,F,G,H]):
                    raise Exception("Invalid segment endpoints in Eqratio")
                ab = _segment_length(A, B)
                cd = _segment_length(C, D)
                ef = _segment_length(E, F)
                gh = _segment_length(G, H)
                ratio1 = (ab / cd) if cd != 0 else float('inf')
                ratio2 = (ef / gh) if gh != 0 else float('inf')
                ok = (cd != 0 and gh != 0 and math.isclose(ratio1, ratio2, rel_tol=1e-6, abs_tol=1e-6))
                if not ok:
                    sname = lambda P,Q: f"{P.name}{Q.name}"
                    lines.append(
                        f"false eqratio: {pred} | {sname(A,B)}={ab:.6f}, {sname(C,D)}={cd:.6f}, {sname(E,F)}={ef:.6f}, {sname(G,H)}={gh:.6f} | ratios: {sname(A,B)}/{sname(C,D)}={ratio1:.6f}, {sname(E,F)}/{sname(G,H)}={ratio2:.6f}"
                    )
            except Exception:
                lines.append(f"false eqratio: {pred}")
            continue

        # Similar triangles and congruent triangles
        # Congruent triangles (ordered only)
        if isinstance(pred, predicate.Contri1):
            t1 = pred.triangle1
            t2 = pred.triangle2
            if not _verify_congruent(t1, t2, ordered=True):
                lines.append(f"false contri1: {pred}")
            continue

        if isinstance(pred, predicate.Contri1):
            t1 = pred.triangle1
            t2 = pred.triangle2
            if not _verify_congruent(t1, t2, ordered=True):
                lines.append(f"false contri1: {pred}")
            continue

        # if isinstance(pred, predicate.Contri2):
        #     t1 = pred.triangle1
        #     t2 = pred.triangle2
        #     if not _verify_congruent(t1, t2, ordered=False):
        #         lines.append(f"false contri2: {pred}")
        #     continue

        # Para: parallel lines AB ∥ CD
        if isinstance(pred, predicate.Para) and len(pts) == 4:
            A, B, C, D = pts
            ok = _parallel(A, B, C, D)
            if not ok:
                lines.append(f"false para: {pred}")
            continue

        # Perp: 4 points P1,P2,P3,P4 representing two rays/segments
        if 'perp' in cls_name and len(pts) == 4:
            P1, P2, P3, P4 = pts
            # angle between vectors (P2-P1) and (P4-P3)
            A1 = _angle_from_points(P2, P1, P3)  # angle at P1 between P2 and P3
            # build vectors
            v1 = np.array([P2.x - P1.x, P2.y - P1.y], dtype=float)
            v2 = np.array([P4.x - P3.x, P4.y - P3.y], dtype=float)
            n1 = np.linalg.norm(v1)
            n2 = np.linalg.norm(v2)
            equal = False
            if n1 > 0 and n2 > 0:
                det = v1[0] * v2[1] - v1[1] * v2[0]
                dot = float(np.dot(v1, v2))
                ang = abs(math.atan2(det, dot))
                # Check near 90° modulo π
                near = min(abs(ang - math.pi/2), abs((math.pi - ang) - math.pi/2))
                equal = near <= tol
            if not equal:
                lines.append(f"perp {P1.name}{P2.name} ⟂ {P3.name}{P4.name} -> is_perp={equal}")
            count += 1
            continue

        # Cong: 4 points A,B,C,D for segments AB and CD
        if 'cong' in cls_name and len(pts) == 4:
            A, B, C, D = pts
            l1 = _segment_length(A, B)
            l2 = _segment_length(C, D)
            equal = abs(l1 - l2) <= (1e-6 if tol < 1e-6 else tol)
            if not equal:
                # Recursive collection of ALL premises trees with sources and rule, not just bad ones
                def collect_premise_tree(p, visited):
                    key = (type(p).__name__, str(p))
                    if key in visited:
                        return { 'pred': p, 'sources': {'cycle'}, 'rule': None, 'ok': None, 'children': [] }
                    visited.add(key)
                    deps_list = []
                    try:
                        deps_list.extend(reason_index_obj.get(p, []) or [])
                    except Exception:
                        pass
                    try:
                        deps_list.extend(reason_index_str.get(str(p), []) or [])
                    except Exception:
                        pass
                    sources = set()
                    rule = None
                    children = []
                    # Compute ok status for this premise
                    try:
                        ok_status = _predicate_ok(p, tol)
                    except Exception:
                        ok_status = False
                    for deps in deps_list:
                        src = deps.get('type') or deps.get('source') or 'unknown'
                        sources.add(str(src))
                        # capture rule when source is DD
                        if rule is None and (str(src).lower().startswith('dd') or str(src).upper() == 'DD'):
                            rule = (deps.get('rule') or deps.get('rule_name') or deps.get('dd_rule') or deps.get('applied_rule'))
                        for sub in deps.get('premises', []) or []:
                            try:
                                children.append(collect_premise_tree(sub, visited))
                            except Exception:
                                children.append({ 'pred': sub, 'sources': {'unknown'}, 'rule': None, 'ok': None, 'children': [] })
                    return { 'pred': p, 'sources': sources if sources else {'unknown'}, 'rule': rule, 'ok': ok_status, 'children': children }

                # Build tree for the false congruence's premises (eqangle, eqratio, simtri, etc.)
                top_sources = set()
                premise_trees = []
                top_deps = []
                try:
                    top_deps.extend(reason_index_obj.get(pred, []) or [])
                except Exception:
                    pass
                try:
                    top_deps.extend(reason_index_str.get(str(pred), []) or [])
                except Exception:
                    pass
                for deps in top_deps:
                    src = deps.get('type') or deps.get('source') or 'unknown'
                    top_sources.add(str(src))
                    for p in deps.get('premises', []) or []:
                        try:
                            premise_trees.append(collect_premise_tree(p, set()))
                        except Exception:
                            premise_trees.append({ 'pred': p, 'sources': {'unknown'}, 'rule': None, 'ok': None, 'children': [] })

                def _format_eqratio_detail(p):
                    try:
                        obj = p
                        if not isinstance(obj, predicate.Eqratio):
                            obj = pred_index_by_str.get(str(p))
                        if obj is None or not isinstance(obj, predicate.Eqratio):
                            return ""
                        segs = obj.ratios
                        r1, r2 = segs[0]
                        r3, r4 = segs[1]
                        A,B = _seg_endpoints(r1)
                        C,D = _seg_endpoints(r2)
                        E,F = _seg_endpoints(r3)
                        G,H = _seg_endpoints(r4)
                        if any(p is None for p in [A,B,C,D,E,F,G,H]):
                            return ""
                        ab = _segment_length(A, B)
                        cd = _segment_length(C, D)
                        ef = _segment_length(E, F)
                        gh = _segment_length(G, H)
                        ratio1 = (ab / cd) if cd != 0 else float('inf')
                        ratio2 = (ef / gh) if gh != 0 else float('inf')
                        sname = lambda P,Q: f"{P.name}{Q.name}"
                        return (f" {sname(A,B)}={ab:.6f}, {sname(C,D)}={cd:.6f}, "
                            f"{sname(E,F)}={ef:.6f}, {sname(G,H)}={gh:.6f} | "
                            f"ratios: {sname(A,B)}/{sname(C,D)}={ratio1:.6f}, "
                            f"{sname(E,F)}/{sname(G,H)}={ratio2:.6f}")
                    except Exception:
                        return ""

                def format_tree(nodes, indent="  "):
                    out = []
                    for n in nodes:
                        p = n.get('pred')
                        srcs = '|'.join(sorted(n.get('sources', {'unknown'})))
                        rule = n.get('rule') or 'none'
                        ok = n.get('ok')
                        status = 'ok' if ok else ('bad' if ok is not None else 'unknown')
                        extra = ""
                        try:
                            if isinstance(p, predicate.Eqratio):
                                extra = _format_eqratio_detail(p)
                        except Exception:
                            extra = ""
                        out.append(f"{indent}- {str(p)} [{status}] | source: {srcs} | rule: {rule}{extra}")
                        kids = n.get('children', [])
                        if kids:
                            out.extend(format_tree(kids, indent + "  "))
                    return out

                src_str = '|'.join(sorted(top_sources)) if top_sources else 'unknown'
                header = f"false cong: {pred} | source: {src_str} | premises:"
                lines.append(header)
                if premise_trees:
                    lines.extend(format_tree(premise_trees, indent="  "))
                else:
                    lines.append("  - none")
            count += 1
            continue

        # Midp: 3 points M,A,B (midpoint of AB)
        if 'midp' in cls_name and len(pts) == 3:
            M, A, B = pts
            # Check distances equal and collinear
            d1 = _segment_length(M, A)
            d2 = _segment_length(M, B)
            equal = abs(d1 - d2) <= (1e-6 if tol < 1e-6 else tol)
            is_col = _collinear(A, B, M)
            ok = equal and is_col
            if not ok:
                lines.append(f"midp {M.name} of {A.name}{B.name} -> is_midpoint={ok} [AM={d1:.6f}, MB={d2:.6f}, collinear={is_col}]")
            count += 1
            continue

        # Sameclock: orientation matching of two triples
        if isinstance(pred, predicate.Sameclock) and len(pts) == 6:
            a1, o1, b1, a2, o2, b2 = pts
            ok = _sameclock(a1, o1, b1, a2, o2, b2)
            if not ok:
                lines.append(f"false sameclock: {pred}")
            continue

        # Cyclic: four points lie on a circle
        if isinstance(pred, predicate.Cyclic) and len(pts) == 4:
            P1, P2, P3, P4 = pts
            center, R = _fit_circle_through_three(P1, P2, P3)
            ok = False
            if center is not None and R is not None:
                ok = _on_circle(center, R, P4)
            if not ok:
                lines.append(f"false cyclic: {pred}")
            continue

        # Circle: center plus three points on the circle
        if isinstance(pred, predicate.Circle) and len(pts) == 4:
            C, P1, P2, P3 = pts
            R = _segment_length(C, P1)
            ok = (_on_circle((C.x, C.y), R, P2) and _on_circle((C.x, C.y), R, P3))
            if not ok:
                lines.append(f"false circle: {pred}")
            continue

    # Use UTF-8 to allow symbols like the perpendicular sign without codec errors
    with open(outfile_path, 'w', encoding='utf-8', errors='replace') as f:
        f.write('\n'.join(lines))

    print(f"Eqangle verification written to {outfile_path} ({len(lines)} counterexamples out of {count} items)")
    return count


def _triangle_side_lengths(t):
    a, b, c = t
    return (
        _segment_length(a, b),
        _segment_length(b, c),
        _segment_length(c, a),
    )


def _verify_congruent(t1, t2, ordered: bool, tol: float = 1e-6) -> bool:
    # Angle checks
    A1 = _angle_from_points(t1[0], t1[1], t1[2])
    B1 = _angle_from_points(t1[1], t1[2], t1[0])
    C1 = _angle_from_points(t1[2], t1[0], t1[1])
    A2 = _angle_from_points(t2[0], t2[1], t2[2])
    B2 = _angle_from_points(t2[1], t2[2], t2[0])
    C2 = _angle_from_points(t2[2], t2[0], t2[1])

    if ordered:
        angles_ok = _equal_mod_pi(A1, A2, tol) and _equal_mod_pi(B1, B2, tol) and _equal_mod_pi(C1, C2, tol)
    else:
        angles_ok = _equal_mod_pi(A1, C2, tol) and _equal_mod_pi(B1, B2, tol) and _equal_mod_pi(C1, A2, tol)

    if not angles_ok:
        return False

    s1 = _triangle_side_lengths(t1)
    s2 = _triangle_side_lengths(t2)
    return (math.isclose(s1[0], s2[0], rel_tol=tol, abs_tol=tol)
            and math.isclose(s1[1], s2[1], rel_tol=tol, abs_tol=tol)
            and math.isclose(s1[2], s2[2], rel_tol=tol, abs_tol=tol))


if __name__ == "__main__":
    # Minimal CLI: try to reuse the last problem by running Solve-like flow
    # This stub demonstrates usage; adapt to your own runner or call verify_eqangles(problem).
    try:
        from Reading_in_Geogebra_File import read_in
        from Problem import GeometricProblem as GP
        from ddar import DDARSystem
        import os

        # Example: use the same test file Solve.py uses (adjust if needed)
        gg_file = os.path.join('Geogebra Files', 'test2.ggb')
        rel_file = os.path.join('test_files', 'problem5_2.txt')

        problem = GP(gg_file, rel_file)
        solver = DDARSystem()
        solver.solve_problem(problem)
        verify_eqangles(problem)
    except Exception as e:
        print(f"verify_eqangles standalone run failed: {e}")