# -*- coding: utf-8 -*-
"""
Created on Tue Apr 28 10:18:07 2026

@author: Bijan Mohammadi

GNU General Public License v3.0

Copyright (C) 2026 Bijan Mohammadi

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

"""

from sympy import *

import json
import os
import sys
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import dblquad
from itertools import combinations
from scipy.integrate import quad
#import random
import warnings
warnings.filterwarnings('ignore')

# n_cpu = os.cpu_count()
# print("n_cpu=",n_cpu)

output_dir = Path('./output')
output_dir.mkdir(exist_ok=True)

verbose=1

SEED = 123
rng = np.random.default_rng(SEED)
GLOBAL_RNG = rng

nfl=20  #nb float digits
nflg="."+str(nfl)+"g"

#%%
def fmt_point(P):
    return [
        format(float(P[0]), nflg),
        format(float(P[1]), nflg)
    ]

def load_analytical_solution(case, k):
    filename = f"Delzant_solution_case_{case}_k"+str(k)+"_best.json"

    if not os.path.exists(output_dir / filename):
        raise FileNotFoundError(
            f"Analytical solution file not found: {filename}"
        )

    with open(output_dir / filename, "r") as f:
        data = json.load(f)

    k = data["k"]

    l = np.array([
        Float(data["l"]["l1"], nfl),
        Float(data["l"]["l2"], nfl),
        Float(data["l"]["l3"], nfl),
        Float(data["l"]["l4"], nfl),
        Float(data["l"]["l5"], nfl)
    ], dtype=float)

    alpha = [Float(a, nfl) for a in data["alpha"]]    

    min_L = float(data["min_L"])
       
    e1_analytical = [
        np.array(P, dtype=float)
        for P in data["geometry"]["e1"]
    ]
    
    e2_analytical = [
        np.array(P, dtype=float)
        for P in data["geometry"]["e2"]
    ]
    
    p1_analytical = np.array(data["geometry"]["p1"], dtype=float)
    p2_analytical = np.array(data["geometry"]["p2"], dtype=float)

    intersections_analytical = [
        np.array(P, dtype=float)
        for P in data["geometry"]["intersections"]
    ]
    
    return k, l, alpha, min_L, e1_analytical, e2_analytical, p1_analytical, p2_analytical, intersections_analytical

#%%
   
plot_poly_and_stop=False
plot_allpoly_and_stop=False

#%%

def octogone_donaldson_delzant_fct(analytical, case, kmin, kmax):

    #analytical=True
    
    if analytical:
        analytical_k_values, analytical_l_values, analytical_alpha_values, analytical_Min_L,\
        e1_analytical, e2_analytical, p1_analytical, p2_analytical, intersections_analytical\
        = load_analytical_solution(case, kmin)
    
        if(kmin!=analytical_k_values):
            print("NOT SAME k analytical")
            sys.exit(0)
        
        kmin=analytical_k_values
        kmax=kmin+1
    
    #search domain min/max for x=[l0,...,l5]
    #l6 is set to 1
    xmin = 1.e-5
    xmax = 300
     
    mc_iter_max = 100
    
    step_size = 1
    num_tests_lines=100
       
    # Numerical integration
    INT_EPS = 1e-14
       
    # Geometry 
    GEOM_EPS = 1e-10
    
    STOP_TOL=-INT_EPS
        
    if(analytical):
        mc_iter_max = 1
        
    #print('INT_EPS=',INT_EPS,'GEOM_EPS=',GEOM_EPS,'STOP_TOL=',STOP_TOL)
        
    #DEFINE PARAMETERS IN MAIN 
    
    #%%
    
    #define different classes
        
    def build_vertices_1(l1, l2, l3, l4, l5, l6, k):
       
        v0 = (0, 0)
       
        v1 = (
            l6 - (k+4)*l1 - (k+3)*l2 - (k+2)*l3 - (k+1)*l4 - k*l5,
            0
        )
    
        v2 = (
            l6 - (k+3)*l2 - (k+2)*l3 - (k+1)*l4 - k*l5,
            l1
        )
    
        v3 = (
            l6 - (k+2)*l3 - (k+1)*l4 - k*l5,
            l1 + l2
        )
    
        v4 = (
            l6 - (k+1)*l4 - k*l5,
            l1 + l2 + l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + l2 + l3 + l4
        )
    
        v6 = (
            l6,
            l1 + l2 + l3 + l4 + l5
        )
    
        v7 = (
            0,
            l1 + l2 + l3 + l4 + l5
        )
    
        # Return as numpy array
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
        
    def build_vertices_2(l1, l2, l3, l4, l5, l6, k):
        v0 = (0, 0)
    
        v1 = (
            l6 - (k+3)*l1 - (k+2)*l2 - (k+1)*l3 - (2*k+1)*l4 - k*l5,
            0
        )
    
        v2 = (
            l6 - (k+2)*l2 - (k+1)*l3 - (2*k+1)*l4 - k*l5,
            l1
        )
    
        v3 = (
            l6 - (k+1)*l3 - (2*k+1)*l4 - k*l5,
            l1 + l2
        )
    
        v4 = (
            l6 - (2*k+1)*l4 - k*l5,
            l1 + l2 + l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + l2 + l3 + 2*l4
        )
    
        v6 = (
            l6,
            l1 + l2 + l3 + 2*l4 + l5
        )
    
        v7 = (
            0,
            l1 + l2 + l3 + 2*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])    
    
    def build_vertices_3(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - (k+3)*l1 - (k+2)*l2 - (2*k+3)*l3 - (k+1)*l4 - k*l5,
            0
        )
    
        v2 = (
            l6 - (k+2)*l2 - (2*k+3)*l3 - (k+1)*l4 - k*l5,
            l1
        )
    
        v3 = (
            l6 - (2*k+3)*l3 - (k+1)*l4 - k*l5,
            l1 + l2
        )
    
        v4 = (
            l6 - (k+1)*l4 - k*l5,
            l1 + l2 + 2*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + l2 + 2*l3 + l4
        )
    
        v6 = (
            l6,
            l1 + l2 + 2*l3 + l4 + l5
        )
    
        v7 = (
            0,
            l1 + l2 + 2*l3 + l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    
    def build_vertices_4(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (2*k+1)*l4 - (3*k+2)*l3 - (k+1)*l2 - (k+2)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (2*k+1)*l4 - (3*k+2)*l3 - (k+1)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (2*k+1)*l4 - (3*k+2)*l3,
            l1 + l2
        )
    
        v4 = (
            l6 - k*l5 - (2*k+1)*l4,
            l1 + l2 + 3*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + l2 + 3*l3 + 2*l4
        )
    
        v6 = (
            l6,
            l1 + l2 + 3*l3 + 2*l4 + l5
        )
    
        v7 = (
            0,
            l1 + l2 + 3*l3 + 2*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])    
    
    def build_vertices_5(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3 - (k+1)*l2 - (k+2)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3 - (k+1)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3,
            l1 + l2
        )
    
        v4 = (
            l6 - k*l5 - (3*k+1)*l4,
            l1 + l2 + 2*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + l2 + 2*l3 + 3*l4
        )
    
        v6 = (
            l6,
            l1 + l2 + 2*l3 + 3*l4 + l5
        )
    
        v7 = (
            0,
            l1 + l2 + 2*l3 + 3*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_6(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (2*k+1)*l4 - (k+1)*l3 - (2*k+3)*l2 - (k+2)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (2*k+1)*l4 - (k+1)*l3 - (2*k+3)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (2*k+1)*l4 - (k+1)*l3,
            l1 + 2*l2
        )
    
        v4 = (
            l6 - k*l5 - (2*k+1)*l4,
            l1 + 2*l2 + l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + 2*l2 + l3 + 2*l4
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + l3 + 2*l4 + l5
        )
    
        v7 = (
            0,
            l1 + 2*l2 + l3 + 2*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_7(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (4*k+1)*l4 - (3*k+1)*l3 - (2*k+1)*l2 - (k+1)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (4*k+1)*l4 - (3*k+1)*l3 - (2*k+1)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (4*k+1)*l4 - (3*k+1)*l3,
            l1 + 2*l2
        )
    
        v4 = (
            l6 - k*l5 - (4*k+1)*l4,
            l1 + 2*l2 + 3*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + 2*l2 + 3*l3 + 4*l4
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + 3*l3 + 4*l4 + l5
        )
    
        v7 = (
            0,
            l1 + 2*l2 + 3*l3 + 4*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_8(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (3*k+1)*l4 - (5*k+2)*l3 - (2*k+1)*l2 - (k+1)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (3*k+1)*l4 - (5*k+2)*l3 - (2*k+1)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (3*k+1)*l4 - (5*k+2)*l3,
            l1 + 2*l2
        )
    
        v4 = (
            l6 - k*l5 - (3*k+1)*l4,
            l1 + 2*l2 + 5*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + 2*l2 + 5*l3 + 3*l4
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + 5*l3 + 3*l4 + l5
        )
    
        v7 = (
            0,
            l1 + 2*l2 + 5*l3 + 3*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
        
    def build_vertices_9(l1, l2, l3, l4, l5, l6, k):
        """
        Build vertices for the given staircase polytope (heavy coefficients case).
    
        Parameters:
            l1,...,l5 : positive parameters
            l6 : typically 1
            k  : parameter
    
        Returns:
            numpy array (8,2)
        """
    
        v0 = (0, 0)
    
        v1 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3 - (3*k+2)*l2 - (k+1)*l1,
            0
        )
    
        v2 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3 - (3*k+2)*l2,
            l1
        )
    
        v3 = (
            l6 - k*l5 - (3*k+1)*l4 - (2*k+1)*l3,
            l1 + 3*l2
        )
    
        v4 = (
            l6 - k*l5 - (3*k+1)*l4,
            l1 + 3*l2 + 2*l3
        )
    
        v5 = (
            l6 - k*l5,
            l1 + 3*l2 + 2*l3 + 3*l4
        )
    
        v6 = (
            l6,
            l1 + 3*l2 + 2*l3 + 3*l4 + l5
        )
    
        v7 = (
            0,
            l1 + 3*l2 + 2*l3 + 3*l4 + l5
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_10(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 + l5 - k*l4 - (k+1)*l3 - (k+2)*l2 - (k+3)*l1,
            0
        )
    
        v2 = (
            l6 + l5 - k*l4 - (k+1)*l3 - (k+2)*l2,
            l1
        )
    
        v3 = (
            l6 + l5 - k*l4 - (k+1)*l3,
            l1 + l2
        )
    
        v4 = (
            l6 + l5 - k*l4,
            l1 + l2 + l3
        )
    
        v5 = (
            l6 + l5,
            l1 + l2 + l3 + l4
        )
    
        v6 = (
            l6,
            l1 + l2 + l3 + l4
        )
    
        v7 = (
            0,
            l1 + l2 + l3 + l4 - l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_11(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 + l5 - k*l4 - (2*k+1)*l3 - (k+1)*l2 - (k+2)*l1,
            0
        )
    
        v2 = (
            l6 + l5 - k*l4 - (2*k+1)*l3 - (k+1)*l2,
            l1
        )
    
        v3 = (
            l6 + l5 - k*l4 - (2*k+1)*l3,
            l1 + l2
        )
    
        v4 = (
            l6 + l5 - k*l4,
            l1 + l2 + 2*l3
        )
    
        v5 = (
            l6 + l5,
            l1 + l2 + 2*l3 + l4
        )
    
        v6 = (
            l6,
            l1 + l2 + 2*l3 + l4
        )
    
        v7 = (
            0,
            l1 + l2 + 2*l3 + l4 - l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])    
    
    def build_vertices_12(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 + l5 - k*l4 - (3*k+1)*l3 - (2*k+1)*l2 - (k+1)*l1,
            0
        )
    
        v2 = (
            l6 + l5 - k*l4 - (3*k+1)*l3 - (2*k+1)*l2,
            l1
        )
    
        v3 = (
            l6 + l5 - k*l4 - (3*k+1)*l3,
            l1 + 2*l2
        )
    
        v4 = (
            l6 + l5 - k*l4,
            l1 + 2*l2 + 3*l3
        )
    
        v5 = (
            l6 + l5,
            l1 + 2*l2 + 3*l3 + l4
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + 3*l3 + l4
        )
    
        v7 = (
            0,
            l1 + 2*l2 + 3*l3 + l4 - l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])    
    
    def build_vertices_13(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 + l5 + l4 - k*l3 - (k+1)*l2 - (k+2)*l1,
            0
        )
    
        v2 = (
            l6 + l5 + l4 - k*l3 - (k+1)*l2,
            l1
        )
    
        v3 = (
            l6 + l5 + l4 - k*l3,
            l1 + l2
        )
    
        v4 = (
            l6 + l5 + l4,
            l1 + l2 + l3
        )
    
        v5 = (
            l6 + l5,
            l1 + l2 + l3
        )
    
        v6 = (
            l6,
            l1 + l2 + l3 - l5
        )
    
        v7 = (
            0,
            l1 + l2 + l3 - l5 - l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
    
    def build_vertices_14(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 + 2*l5 + l4 - k*l3 - (2*k+1)*l2 - (k+1)*l1,
            0
        )
    
        v2 = (
            l6 + 2*l5 + l4 - k*l3 - (2*k+1)*l2,
            l1
        )
    
        v3 = (
            l6 + 2*l5 + l4 - k*l3,
            l1 + 2*l2
        )
    
        v4 = (
            l6 + 2*l5 + l4,
            l1 + 2*l2 + l3
        )
    
        v5 = (
            l6 + 2*l5,
            l1 + 2*l2 + l3
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + l3 - l5
        )
    
        v7 = (
            0,
            l1 + 2*l2 + l3 - l5 - l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])
        
    
    def build_vertices_15(l1, l2, l3, l4, l5, l6, k):
    
        v0 = (0, 0)
    
        v1 = (
            l6 - (k+1)*l1 - (2*k+1)*l2 - k*l3 + l4 + l5,
            0
        )
    
        v2 = (
            l6 - (2*k+1)*l2 - k*l3 + l4 + l5,
            l1
        )
    
        v3 = (
            l6 - k*l3 + l4 + l5,
            l1 + 2*l2
        )
    
        v4 = (
            l6 + l4 + l5,
            l1 + 2*l2 + l3
        )
    
        v5 = (
            l6 + l5,
            l1 + 2*l2 + l3
        )
    
        v6 = (
            l6,
            l1 + 2*l2 + l3 - l5
        )
    
        v7 = (
            0,
            l1 + 2*l2 + l3 - l5 - 2*l6
        )
    
        return np.array([v0, v1, v2, v3, v4, v5, v6, v7])    
    
    #%%
    def get_primitive_dirs(case, k):
        """
        Primitive directions of the 8 edges.
    
        Edge ordering:
          0 : bottom horizontal
          1..5 : staircase edges
          6 : top horizontal (leftward)
          7 : left vertical (downward)
        """
    
        case = str(case)
    
        dirs = {
            "1": [
                (1,0),
                (k+4,1),
                (k+3,1),
                (k+2,1),
                (k+1,1),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "2": [
                (1,0),
                (k+3,1),
                (k+2,1),
                (k+1,1),
                (2*k+1,2),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "3": [
                (1,0),
                (k+3,1),
                (k+2,1),
                (2*k+3,2),
                (k+1,1),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "4": [
                (1,0),
                (k+2,1),
                (k+1,1),
                (3*k+2,3),
                (2*k+1,2),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "5": [
                (1,0),
                (k+2,1),
                (k+1,1),
                (2*k+1,2),
                (3*k+1,3),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "6": [
                (1, 0),
                (k+2, 1),
                (2*k+3, 2),
                (k+1, 1),
                (2*k+1, 2),
                (k, 1),
                (-1, 0),
                (0, -1),
            ],        
    
            "7": [
                (1,0),
                (k+1,1),
                (2*k+1,2),
                (3*k+1,3),
                (4*k+1,4),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "8": [
                (1,0),
                (k+1,1),
                (2*k+1,2),
                (5*k+2,5),
                (3*k+1,3),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "9": [
                (1,0),
                (k+1,1),
                (3*k+2,3),
                (2*k+1,2),
                (3*k+1,3),
                (k,1),
                (-1,0),
                (0,-1),
            ],
    
            "10": [
                (1,0),
                (k+3,1),
                (k+2,1),
                (k+1,1),
                (k,1),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
    
            "11": [
                (1,0),
                (k+2,1),
                (k+1,1),
                (2*k+1,2),
                (k,1),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
    
            "12": [
                (1,0),
                (k+1,1),
                (2*k+1,2),
                (3*k+1,3),
                (k,1),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
    
            "13": [
                (1,0),
                (k+2,1),
                (k+1,1),
                (k,1),
                (-1,0),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
    
            "14": [
                (1,0),
                (k+1,1),
                (2*k+1,2),
                (k,1),
                (-1,0),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
    
            "15": [
                (1,0),
                (k+1,1),
                (2*k+1,2),
                (k,1),
                (-1,0),
                (-1,0),
                (-1,0),
                (0,-1),
            ],
        }
    
        return dirs[case]
    
    #%%
    def is_convex_polygon(points):
    
        n = len(points)
        prev = 0
    
        for i in range(n):
            p0 = points[i]
            p1 = points[(i+1)%n]
            p2 = points[(i+2)%n]
    
            v1 = p1 - p0
            v2 = p2 - p1
    
            cross = v1[0]*v2[1] - v1[1]*v2[0]
    
            if abs(cross) > GEOM_EPS:
                if prev == 0:
                    prev = np.sign(cross)
                else:
                    if np.sign(cross) != prev:
                        return False
    
        return True
    
    #%%
    if(case=="1"):
        build_vertices=build_vertices_1
    elif(case=="2"):
        build_vertices=build_vertices_2
    elif(case=="3"):
        build_vertices=build_vertices_3
    elif(case=="4"):
        build_vertices=build_vertices_4
    elif(case=="5"):
        build_vertices=build_vertices_5
    elif(case=="6"):
        build_vertices=build_vertices_6
    elif(case=="7"):
        build_vertices=build_vertices_7
    elif(case=="8"):
        build_vertices=build_vertices_8
    elif(case=="9"):
        build_vertices=build_vertices_9
    elif(case=="10"):
        build_vertices=build_vertices_10
    elif(case=="11"):
        build_vertices=build_vertices_11
    elif(case=="12"):
        build_vertices=build_vertices_12
    elif(case=="13"):
        build_vertices=build_vertices_13
    elif(case=="14"):
        build_vertices=build_vertices_14
    elif(case=="15"):
        build_vertices=build_vertices_15
    else:
        print("no vertices fct")
        sys.exit(0)
    
    #%%
    
    def func_L(x, k):
        
        l1=x[0]
        l2=x[1]
        l3=x[2]
        l4=x[3]
        l5=x[4]
        
        semiglobal_min_functional=np.inf
        semiglobal_geom = None
    
        l6=1
        verts=build_vertices(l1, l2, l3, l4, l5, l6, k)
        
        if is_convex_polygon(verts) == False:
            return 10000, None, None
        
        #test polygone shape
        if(plot_poly_and_stop):
            V = np.vstack([verts, verts[0]])
            print(V)
            plt.plot(V[:,0], V[:,1], 'o-')
            plt.axis('equal')
            plt.show()
            sys.exit(0)
       
        n = len(verts)
        edges = [(verts[i], verts[(i + 1) % n]) for i in range(n)]
        
        if(False): #geometry test Delzant polygone
            for a,b in edges:
                dx = b[0]-a[0]
                dy = b[1]-a[1]
                print(dx,dy)
                
            for i,(a,b) in enumerate(edges):        
                dx = b[0]-a[0]
                dy = b[1]-a[1]        
                L = np.hypot(dx,dy)
                print(i, dx, dy, L/abs(dx))   
            sys.exit(0)
    
    #%%
    #ALL NECESSARY FUNCTIONS
    
        primitive_dirs = get_primitive_dirs(case, k)
        
        
        def edge_sigma_weight(edge_index):
            vx, vy = primitive_dirs[edge_index]
        
            a, b = edges[edge_index]
        
            dx = b[0] - a[0]
            dy = b[1] - a[1]
        
            if abs(vx) > 0:
                m = dx / vx
            else:
                m = dy / vy
        
            return abs(m)
            
        def poly_integral_of_affine(c0, c1, c2):
            xmin, xmax = verts[:, 0].min(), verts[:, 0].max()
            def integrand(x2, x1):
                return c0 + c1 * x1 + c2 * x2
            return dblquad(integrand, xmin, xmax, lambda x1: bounds_for_x(x1)[0], 
                           lambda x1: bounds_for_x(x1)[1], epsabs=INT_EPS, epsrel=INT_EPS)[0]
        
        def poly_moment(i, j):
            xmin, xmax = verts[:, 0].min(), verts[:, 0].max()
        
            def integrand(x2, x1):
                vals = [1.0, x1, x2]
                return vals[i] * vals[j]
        
            return dblquad(
                integrand,
                xmin,
                xmax,
                lambda x1: bounds_for_x(x1)[0],
                lambda x1: bounds_for_x(x1)[1],
                epsabs=INT_EPS,
                epsrel=INT_EPS,
            )[0]
        
        def boundary_moment_euclidean(i):
            total = 0.0
            for a, b in edges:
                L = np.linalg.norm(b - a)
        
                def integrand(t):
                    x = a[0] + t * (b[0] - a[0])
                    y = a[1] + t * (b[1] - a[1])
                    vals = [1.0, x, y]
                    return vals[i]
            
                total += L * quad(integrand, 0.0, 1.0, epsabs=INT_EPS, epsrel=INT_EPS)[0]
            return total
        
        def boundary_moment_delzant(i):
            total = 0.0
        
            for edge_index, (a, b) in enumerate(edges):
        
                sigma = edge_sigma_weight(edge_index)
        
                def integrand(t):
                    x = a[0] + t * (b[0] - a[0])
                    y = a[1] + t * (b[1] - a[1])
        
                    if i == 0:
                        return 1.0
                    elif i == 1:
                        return x
                    else:
                        return y
        
                val, _ = quad(
                    integrand,
                    0.0,
                    1.0,
                    epsabs=INT_EPS,
                    epsrel=INT_EPS,
                )
        
                total += sigma * val
        
            return total    
        
        boundary_moment = boundary_moment_delzant
        
        def solve_alpha_from_polygon_moments():
            M = np.array([
                [poly_moment(0, 0), poly_moment(0, 1), poly_moment(0, 2)],
                [poly_moment(1, 0), poly_moment(1, 1), poly_moment(1, 2)],
                [poly_moment(2, 0), poly_moment(2, 1), poly_moment(2, 2)],
            ], dtype=float)
        
            b = np.array([
                boundary_moment(0),
                boundary_moment(1),
                boundary_moment(2),
            ], dtype=float)
        
            alpha = np.linalg.solve(M, b)
            info = {
                "M": M,
                "b": b,
                "residual": M @ alpha - b,
            }
            # print("cond=",np.linalg.cond(info["M"]))
            # sys.exit(0)
            
            return alpha, info
        
        def check_normalization(alpha, tol=GEOM_EPS):
            M = np.array([
                [poly_moment(0, 0), poly_moment(0, 1), poly_moment(0, 2)],
                [poly_moment(1, 0), poly_moment(1, 1), poly_moment(1, 2)],
                [poly_moment(2, 0), poly_moment(2, 1), poly_moment(2, 2)],
            ], dtype=float)
        
            b = np.array([
                boundary_moment(0),
                boundary_moment(1),
                boundary_moment(2),
            ], dtype=float)
        
            residual = M @ alpha - b
            ok = np.all(np.abs(residual) < tol)
            return ok, residual.tolist()
        
        def cross2(a, b):
            return a[0] * b[1] - a[1] * b[0]
        
        
        def polygon_area(poly):
            x = poly[:, 0]
            y = poly[:, 1]
            return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
        
        def bounds_for_x(x):
            #n is nbre of vertices : 8
            ys = []
            for i in range(n):
                a = verts[i]
                b = verts[(i + 1) % n]
                if abs(a[0] - b[0]) < GEOM_EPS:
                    if abs(x - a[0]) < GEOM_EPS:
                        ys.extend([a[1], b[1]])
                    continue
                if x < min(a[0], b[0]) - GEOM_EPS or x > max(a[0], b[0]) + GEOM_EPS:
                    continue
                t = (x - a[0]) / (b[0] - a[0])
                if -GEOM_EPS <= t <= 1 + GEOM_EPS:
                    ys.append(a[1] + t * (b[1] - a[1]))
            if len(ys) < 2:
                return None
            return min(ys), max(ys)
        
        ALPHA, info = solve_alpha_from_polygon_moments()
        
        if analytical:
            ALPHA = analytical_alpha_values
        
        ok, vals = check_normalization(ALPHA)
        if(verbose>3):
            print("normalization_ok:", ok)
            print("affine residuals:", vals)
        
        
        def line_from_points(p1, p2):
            dx = p2[0] - p1[0]
            dy = p2[1] - p1[1]
            a = dy
            b = -dx
            c = dx * p1[1] - dy * p1[0]
            return a, b, c
        
        
        def point_on_edge(k, rng):
            a, b = edges[k]
            t = rng.random()
            return a + t * (b - a)
        
        def segment_intersection_interior(p1, p2, q1, q2, eps=GEOM_EPS):
            r = p2 - p1
            s = q2 - q1
            den = cross2(r, s)
            if abs(den) < eps:
                return False
            t = cross2(q1 - p1, s) / den
            u = cross2(q1 - p1, r) / den
            return (eps < t < 1 - eps) and (eps < u < 1 - eps)
        
        
        def hit_edges(p1, p2):
            hits = []
            for k, (q1, q2) in enumerate(edges):
                if segment_intersection_interior(p1, p2, q1, q2):
                    hits.append(k)
            return sorted(set(hits))
        
        
        def edge_boundary_integral_euclidean(u):
            if(False):
                for i in range(len(edges)):
                    w = edge_sigma_weight(i)        
                    a,b = edges[i]
                    L = np.linalg.norm(b-a)
                    print(i, "sigma=", w, "euclid=", L, "ratio=", L/w)            
                sys.exit(0)
    
            total = 0.0
            for a, b in edges:
                L = np.linalg.norm(b - a)
                def integrand(t):
                    x = a[0] + t * (b[0] - a[0])
                    y = a[1] + t * (b[1] - a[1])
                    return u(x, y)
                total += L * quad(integrand, 0.0, 1.0, epsabs=INT_EPS, epsrel=INT_EPS)[0]
            return total
        
        def edge_boundary_integral_delzant(u):
            total = 0.0
        
            for edge_index, (a, b) in enumerate(edges):
        
                sigma = edge_sigma_weight(edge_index)
        
                def integrand(t):
                    x = a[0] + t * (b[0] - a[0])
                    y = a[1] + t * (b[1] - a[1])
                    return u(x, y)
        
                total += sigma * quad(
                    integrand,
                    0.0,
                    1.0,
                    epsabs=INT_EPS,
                    epsrel=INT_EPS
                )[0]
        
            return total    
        
        edge_boundary_integral = edge_boundary_integral_delzant
        
        def compute_LA(p1, p2, alpha=ALPHA):
            globals()["compute_LA_global"] = compute_LA
            a, b, c = line_from_points(p1, p2)
        
            def u(x1, x2):
                return max(0.0, a * x1 + b * x2 + c)
        
            boundary = edge_boundary_integral(u)
        
            def int_Au(x2, x1):
                return (alpha[0] + alpha[1] * x1 + alpha[2] * x2) * u(x1, x2)
        
            xmin, xmax = verts[:, 0].min(), verts[:, 0].max()
            area_term = dblquad(
                int_Au,
                xmin,
                xmax,
                lambda x1: bounds_for_x(x1)[0],
                lambda x1: bounds_for_x(x1)[1],
                epsabs=INT_EPS,
                epsrel=INT_EPS,
            )[0]
    
            if(False): #check integration
                print("boundary =", boundary)
                print("area_term =", area_term)
                print("LA =", boundary - area_term)
                print("boundary/area =", boundary/area_term)
                sys.exit(0)
        
            return boundary - area_term, boundary, area_term, (a, b, c)
        
        def line_edge_intersections(p1, p2, eps=GEOM_EPS):
            r = p2 - p1
            hits = []
            for k, (q1, q2) in enumerate(edges):
                s = q2 - q1
                den = cross2(r, s)
                if abs(den) < eps:
                    continue
                t = cross2(q1 - p1, s) / den
                u = cross2(q1 - p1, r) / den
                if -eps <= u <= 1 + eps:
                    x = p1 + t * r
                    if min(q1[0], q2[0]) - eps <= x[0] <= max(q1[0], q2[0]) + eps and \
                       min(q1[1], q2[1]) - eps <= x[1] <= max(q1[1], q2[1]) + eps:
                        if eps < u < 1 - eps:
                            hits.append((k, u, x))
            return hits
        
        def is_degenerate(p1, p2, eps=GEOM_EPS):
            hits = line_edge_intersections(p1, p2, eps=eps)
            if len(hits) != 2:
                return True
            return hits[0][0] == hits[1][0]
           
        
        if(plot_allpoly_and_stop):
            def plot_all_cases(df):
                m = len(df)
                fig, axes = plt.subplots(4, 7, figsize=(24, 13), constrained_layout=True)
                axes = axes.ravel()
                poly_closed = np.vstack([verts, verts[0]])
                for ax, (_, row) in zip(axes, df.iterrows()):
                    ax.plot(poly_closed[:, 0], poly_closed[:, 1], 'k-', lw=2)
                    ax.fill(verts[:, 0], verts[:, 1], color='lightblue', alpha=0.35)
                    if abs(row['b']) > GEOM_EPS:
                        xs = np.linspace(verts[:, 0].min() - 0.1, verts[:, 0].max() + 0.1, 300)
                        ys = (-row['a'] * xs - row['c']) / row['b']
                        if(row['LA']<0):
                            ax.plot(xs, ys, color='tab:red', lw=1.8)
                        else:
                            ax.plot(xs, ys, color='tab:green', lw=1.8)
                    ax.set_title(f"{int(row['edge_i'])}×{int(row['edge_j'])} LA={row['LA']:.4f}", fontsize=20)
                    ax.set_xlim(verts[:, 0].min() - 0.15, verts[:, 0].max() + 0.15)
                    ax.set_ylim(verts[:, 1].min() - 0.15, verts[:, 1].max() + 0.15)
                    ax.set_aspect('equal', adjustable='box')
                    ax.grid(True, alpha=0.2)
                    ax.tick_params(axis='both', which='both', labelbottom=False, labelleft=False, bottom=False, left=False)
                for ax in axes[m:]:
                    ax.axis('off')
                    
                fig.suptitle('PL test families-', fontsize=20)    
                fig.savefig(output_dir / "test_families.png", dpi=200)
                plt.show()
        
            plot_rng = np.random.default_rng(0)   
                  
            rows = []
            for idx, (i, j) in enumerate(combinations(range(n), 2), 1):
                p1 = point_on_edge(i, plot_rng)
                p2 = point_on_edge(j, plot_rng)
                deg = is_degenerate(p1, p2)
                
                if deg:
                    continue
                
                LA, iu, iAu, (a, b, c) = compute_LA(p1, p2)
                            
                rows.append({
                    'case': idx,
                    'edge_i': i,
                    'edge_j': j,
                    'p1x': p1[0], 'p1y': p1[1],
                    'p2x': p2[0], 'p2y': p2[1],
                    'a': a, 'b': b, 'c': c,
                    'int_u': iu,
                    'int_Au': iAu,
                    'LA': LA,
                    'degenerate': deg,
                })
            df = pd.DataFrame(rows)
            plot_all_cases(df)
            print(f'Polygon area: {polygon_area(verts):.6f}')
            print(f'Alpha: {ALPHA.tolist()}') 
            min_LA=np.min(df['LA'])
            print(f'Min L_A: {min_LA}') 
            sys.exit(0)
        
        # ============================================================
        #   GÉOMÉTRIE : SEGMENT – SEGMENT INTERSECTION
        # ============================================================
    
        def segment_intersection(p1, p2, q1, q2, eps=GEOM_EPS):
            p = p1
            r = p2 - p1
            q = q1
            s = q2 - q1
        
            rxs = cross2(r, s)
            if abs(rxs) < eps:
                return None
        
            t = cross2(q - p, s) / rxs
            u = cross2(q - p, r) / rxs
        
            if 0 < t < 1 and 0 < u < 1:
                return p + t*r
        
            return None
        
        # ============================================================
        #   DeTECTION : true intersection with 2 distinct edges 
        # ============================================================
        
        def detect_valid_line(p_pt1, p_pt2, verts):
            """
            detects if line by p_pt1–p_pt2 cuts 2 edges interior of octagon        
            """
            n = len(verts)
            intersections = []
            edge_ids = []
        
            for i in range(n):
                A = verts[i]
                B = verts[(i + 1) % n]
                P = segment_intersection(p_pt1, p_pt2, A, B)
                if P is not None:
                    intersections.append(P)
                    edge_ids.append(i)
        
            if len(intersections) != 2:
                return False, intersections, edge_ids
        
            e1, e2 = edge_ids
            if e1 == e2:
                return False, intersections, edge_ids
        
            return True, intersections, edge_ids
        
        
        # ============================================================
        #   DeTECTION : Area  P₊ = {u > 0}
        # ============================================================
        
        def area_P_plus(a, b, c, pts):
            """
            Approx area P₊ = {u > 0} via sampling
            pts = uniform pts in octagon
            """
            X = pts[:, 0]
            Y = pts[:, 1]
            uvals = a * X + b * Y + c
            mask = uvals > 0
            ratio = np.mean(mask)
        
            xs = verts[:, 0]
            ys = verts[:, 1]
            area = 0.5 * abs(np.dot(xs, np.roll(ys, -1)) - np.dot(ys, np.roll(xs, -1)))
        
            return ratio * area
        
        
        # ============================================================
        #  Degeneracy test
        # ============================================================
        
    
        def is_degenerate_line(p1, p2, verts):
            ok, intersections, edge_ids = detect_valid_line(p1, p2, verts)
        
            if not ok:
                return True, "Invalid intersection", intersections
        
            return False, "Valid", intersections
            
        def point_in_polygon(point, poly):
            x, y = point
            inside = False
            n = len(poly)
            j = n - 1
        
            for i in range(n):
                xi, yi = poly[i]
                xj, yj = poly[j]
        
                intersect = ((yi > y) != (yj > y)) and (
                    x < (xj - xi) * (y - yi) / (yj - yi + GEOM_EPS) + xi
                )
                if intersect:
                    inside = not inside
                j = i
        
            return inside
            
        # ==========================================================
        #   TEST N random lines
        # ==========================================================
        
    
        def test_random_lines(edges, verts, compute_LA, num_tests=100, ALPHA=None):
                   
            #rng = random.Random(123)
            rng = GLOBAL_RNG
                        
            local_min = np.inf
            local_geom = None
        
            for kk in range(num_tests):
    
                # e1, e2 = rng.sample(edges, 2)
                
                # p1 = e1[0] + rng.random()*(e1[1]-e1[0])
                # p2 = e2[0] + rng.random()*(e2[1]-e2[0])
                
                idx = rng.choice(len(edges), size=2, replace=False)
    
                e1 = edges[idx[0]]
                e2 = edges[idx[1]]
    
                p1 = e1[0] + rng.random() * (e1[1] - e1[0])
                p2 = e2[0] + rng.random() * (e2[1] - e2[0])
                         
                # heck degeneracy FIRST
                deg, reason, inters = is_degenerate_line(p1, p2, verts)
        
                if deg:
                    #print("DEG:",deg)
                    continue
        
                # only compute L_A for valid lines
    
                # test of reproducibility            
                # for i in range(10):
                #     LA, *_ = compute_LA(p1, p2)
                #     print("{:.18e}".format(LA))
                
                L_A_val, *_ = compute_LA(p1, p2)
        
                if L_A_val < local_min:
                    local_min = L_A_val
                    local_geom = [e1, e2, p1, p2, deg, reason, inters]
        
                # if(L_A_val<0):
                #     print(f"iter {kk} | L_A={L_A_val:.3e} | deg={deg}")
        
                if( (analytical==False and local_min < STOP_TOL) or 
                    (analytical and abs(local_min-analytical_Min_L)<INT_EPS)
                    ):
                    break
        
            return local_min, local_geom
                
        # Lancer le test
        
        if not analytical:
        
            local_min_functional, local_geom = test_random_lines(
                edges, verts, compute_LA,
                num_tests=num_tests_lines,
                ALPHA=ALPHA
            )
        
        else:
        
            LA, B, A, abc = compute_LA(
                p1_analytical,
                p2_analytical,
                alpha=ALPHA
            )
                
            local_min_functional = LA
        
            local_geom = [
                e1_analytical,           # e1
                e2_analytical,           # e2
                p1_analytical,  # p1
                p2_analytical,  # p2
                False,          # degenerate
                "analytical",   # reason
                intersections_analytical              # intersections
            ]
               
                
        if(local_min_functional<semiglobal_min_functional):
            semiglobal_min_functional=local_min_functional
            semiglobal_geom=local_geom
            
        #a posteriori test for analytical solution
        if analytical and semiglobal_geom is not None:
            p1 = semiglobal_geom[2]
            p2 = semiglobal_geom[3]            
            LA_check, *_ = compute_LA(p1, p2)
            if(False):
                print("check numerical integration")
                print("stored     =", semiglobal_min_functional)
                print("recomputed =", LA_check)
                print("diff       =", semiglobal_min_functional - LA_check)        
            
        return semiglobal_min_functional, semiglobal_geom, ALPHA
    
    #end func_L
    
    #%%
    
    #sampling of l_i, l6=1
    
    def compute_S_1(l, k):
        return (k+4)*l[0] + (k+3)*l[1] + (k+2)*l[2] + (k+1)*l[3] + k*l[4] - 1
    
    def compute_S_2(l, k):
        return (k+3)*l[0] + (k+2)*l[1] + (k+1)*l[2] + (2*k+1)*l[3] + k*l[4] - 1
    
    def compute_S_3(l, k):
        return (k+3)*l[0] + (k+2)*l[1] + (2*k+3)*l[2] + (k+1)*l[3] + k*l[4] - 1
    
    def compute_S_4(l, k):
        return -(1 - k*l[4] - (2*k+1)*l[3] - (3*k+2)*l[2] - (k+1)*l[1] - (k + 2)*l[0])
    
    def compute_S_5(l, k):
    #-(ℓ6 − 𝑘ℓ5 −(3𝑘 + 1)ℓ4 − (2𝑘 + 1)ℓ3 − (𝑘 + 1)ℓ2 − (𝑘 + 2)ℓ1) <0 
        return -(1 - k*l[4] - (3*k+1)*l[3] - (2*k+1)*l[2] - (k+1)*l[1] - (k + 2)*l[0])
    
    def compute_S_6(l, k):
    # ℓ6 − 𝑘ℓ5 −(2𝑘 + 1)ℓ4 − (𝑘 + 1)ℓ3 − (2𝑘 + 3)ℓ2 − (𝑘 + 2)ℓ1 > 0
        return -(1 - k*l[4] - (2*k + 1)*l[3] - (k + 1)*l[2] - (2*k + 3)*l[1] - (k + 2)*l[0])
    
    def compute_S_7(l, k):
    # ℓ6 − 𝑘ℓ5 −(4𝑘 + 1)ℓ4 − (3𝑘 + 1)ℓ3 − (2𝑘 + 1)ℓ2 − (𝑘 + 1)ℓ1 > 0
        return -(1 - k*l[4] - (4*k + 1)*l[3] - (3*k + 1)*l[2] - (2*k + 1)*l[1] - (k + 1)*l[0])
    
    def compute_S_8(l, k):
    #ℓ6 − 𝑘ℓ5 −(3𝑘 + 1)ℓ4 − (5𝑘 + 2)ℓ3 − (2𝑘 + 1)ℓ2 − (𝑘 + 1)ℓ1 > 0
        return -(1 - k*l[4] - (3*k + 1)*l[3] - (5*k + 2)*l[2] - (2*k + 1)*l[1] - (k + 1)*l[0])
    
    def compute_S_9(l, k):
    #ℓ6 − 𝑘ℓ5 − (3𝑘 + 1)ℓ4 − (2𝑘 + 1)ℓ3 − (3𝑘 + 2)ℓ2 − (𝑘 + 1)ℓ1 > 0
        return -(1 - k*l[4] - (3*k + 1)*l[3] - (2*k + 1)*l[2] - (3*k + 2)*l[1] - (k + 1)*l[0])
    
    def compute_S_10(l, k):
    #-(ℓ6 + ℓ5 −𝑘ℓ4 − (𝑘 + 1)ℓ3 − (𝑘 + 2)ℓ2 − (𝑘 + 3)ℓ1) < 0 and -(ℓ1 + ℓ2 + ℓ3 + ℓ4 − ℓ6) < 0
        l1, l2, l3, l4, l5 = l
        C1 = -(1 + l5 - k*l4 - (k+1)*l3 - (k+2)*l2 - (k+3)*l1)
        C2 = -(l1 + l2 + l3 + l4 - 1)
        return C1,C2
    
    def compute_S_11(l, k):
    #-(ℓ6 + ℓ5 − 𝑘ℓ4 − (2𝑘 + 1)ℓ3 − (𝑘 + 1)ℓ2 − (𝑘 + 2)ℓ1) < 0 and -(ℓ1 + ℓ2 + 2ℓ3 + ℓ4 − ℓ6) < 0
        l1, l2, l3, l4, l5 = l
        C1 = -(1 + l5 - k*l4 - (2*k+1)*l3 - (k+1)*l2 - (k+2)*l1)
        C2 = -(l1 + l2 + 2*l3 + l4 - 1)
        return C1,C2
    
    def compute_S_12(l, k):
    #-(ℓ6 + ℓ5 − 𝑘ℓ4 − (3𝑘 + 1)ℓ3 − (2𝑘 + 1)ℓ2 − (𝑘 + 1)ℓ1) < 0 and -(ℓ1 + 2ℓ2 + 3ℓ3 + ℓ4 − ℓ6) < 0
        l1, l2, l3, l4, l5 = l
        C1 = -(1 + l5 - k*l4 - (3*k+1)*l3 - (2*k+1)*l2 - (k+1)*l1)
        C2 = -(l1 + 2*l2 + 3*l3 + l4 - 1)
        return C1,C2
    
    def compute_S_13(l, k):
    #-(ℓ6 + 2ℓ5 +ℓ4 − 𝑘ℓ3 − (𝑘 + 1)ℓ2 − (𝑘 + 2)ℓ1) < 0 and -(ℓ1 + ℓ2 + ℓ3 − ℓ5 − ℓ6) < 0
        l1, l2, l3, l4, l5 = l
        C1 = -(1 + 2*l5 + l4 - k*l3 - (k+1)*l2 - (k+2)*l1)
        C2 = -(l1 + l2 + l3 - l5 - 1)
        return C1,C2
    
    def compute_S_14(l, k):
    #ℓ6 + 2ℓ5 +ℓ4 − 𝑘ℓ3 − (2𝑘 + 1)ℓ2 − (𝑘 + 1)ℓ1 > 0 and ℓ1 + 2ℓ2 + ℓ3 − ℓ5 − ℓ6 > 0
        l1, l2, l3, l4, l5 = l
        C1 = -(1 + 2*l5 + l4 - k*l3 - (2*k+1)*l2 - (k+1)*l1)
        C2 = -(l1 + 2*l2 + l3 - l5 - 1)
        return C1,C2
    
    def compute_S_15(l, k):
    #(𝑘 + 1)ℓ1 + (2𝑘 + 1)ℓ2 +𝑘ℓ3 − ℓ4 − ℓ5 − ℓ6 < 0 and ℓ1 + 2ℓ2 + ℓ3 − ℓ5 − 2ℓ6 > 0
        C1=(k+1)*l[0]+(2*k+1)*l[1]+k*l[2]-l[3]-l[4]-1
        C2=-(l[0]+2*l[1]+l[2]-l[4]-2)
        return C1,C2
    
    if(case=="1"):
        compute_S=compute_S_1
    elif(case=="2"):
        compute_S=compute_S_2
    elif(case=="3"):
        compute_S=compute_S_3
    elif(case=="4"):
        compute_S=compute_S_4
    elif(case=="5"):
        compute_S=compute_S_5
    elif(case=="6"):
        compute_S=compute_S_6
    elif(case=="7"):
        compute_S=compute_S_7
    elif(case=="8"):
        compute_S=compute_S_8
    elif(case=="9"):
        compute_S=compute_S_9
    else:
        compute_S=None #for case>9 2 constraints 
    
    #%%
    
    def transform_1(y, k):
        y = y / (1 + np.sum(y))   # enforce constraint            
        return np.array([
            y[0]/(k+4),
            y[1]/(k+3),
            y[2]/(k+2),
            y[3]/(k+1),
            y[4]/(k if k>0 else 1),
        ])
    def transform_2(y, k):
        y = y / (1 + np.sum(y))   # simplex constraint
        return np.array([
            y[0]/(k+3),
            y[1]/(k+2),
            y[2]/(2*k+3),
            y[3]/(k+1),
            y[4]/k if k != 0 else y[4]   # avoid division by 0
        ])
    def transform_3(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0]/(k+3),
            y[1]/(k+2),
            y[2]/(k+1),
            y[3]/(2*k+1),
            y[4]/k if k != 0 else y[4]
        ])
    def transform_4(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+2),
            y[1] / (k+1),
            y[2] / (3*k + 2),
            y[3] / (2*k + 1),
            y[4] / k if k != 0 else y[4]
        ])
    def transform_5(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+2),
            y[1] / (k+1),
            y[2] / (2*k+1),
            y[3] / (3*k+1),
            y[4] / k if k != 0 else y[4]
        ])
    
    def transform_6(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+2),
            y[1] / (2*k+3),
            y[2] / (k+1),
            y[3] / (2*k+1),
            y[4] / k if k != 0 else y[4]
        ])
    def transform_7(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+1),
            y[1] / (2*k+1),
            y[2] / (3*k+1),
            y[3] / (4*k+1),
            y[4] / k if k != 0 else y[4]
        ])
    def transform_8(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+1),
            y[1] / (2*k+1),
            y[2] / (5*k+2),
            y[3] / (3*k+1),
            y[4] / k if k != 0 else y[4]
        ])
    def transform_9(y, k):
        y = y / (1 + np.sum(y))
        return np.array([
            y[0] / (k+1),
            y[1] / (3*k + 2),
            y[2] / (2*k + 1),
            y[3] / (3*k + 1),
            y[4] / k if k != 0 else y[4]
        ])
    
    if(case=="1"):
        transform=transform_1
    elif(case=="2"):
        transform=transform_2
    elif(case=="3"):
        transform=transform_3
    elif(case=="4"):
        transform=transform_4
    elif(case=="5"):
        transform=transform_5
    elif(case=="6"):
        transform=transform_6
    elif(case=="7"):
        transform=transform_7
    elif(case=="8"):
        transform=transform_8
    elif(case=="9"):
        transform=transform_9    
    elif(int(case)>9):
        transform=None  #for 2 constraints, do not use transform in y
    
    #%%
    def sample_near_boundary_1to9(k):
        while True:
            l = np.exp(rng.standard_normal(5))
            S = compute_S(l, k)
            if S <= GEOM_EPS:
                continue
            target = rng.uniform(0.95, 0.99)
            return (target / S) * l
       
    def sample_near_boundary_10(k, max_trials=1000):
    #ℓ6 + ℓ5 − 𝑘ℓ4 − (𝑘 + 1)ℓ3 − (𝑘 + 2)ℓ2 − (𝑘 + 3)ℓ1 > 0 and ℓ1 + ℓ2 + ℓ3 + ℓ4 − ℓ6 > 0, 
        for _ in range(max_trials):
    
            # --- Step 1: sample l1,l2,l3,l4
            #l1234 = np.exp(rng.standard_normal(4))
            l1234 = rng.uniform(0.1, 2.0, size=4)
    
            l1, l2, l3, l4 = l1234
    
            # enforce C2 > 0
            if l1 + l2 + l3 + l4 <= 1:
                continue
    
            # --- Step 2: compute lower bound for l5
            A = k*l4 + (k+1)*l3 + (k+2)*l2 + (k+3)*l1 - 1
    
            # choose l5
            if A <= 0:
                l5 = np.exp(rng.standard_normal()) * 0.2
            else:
                l5 = A + rng.uniform(0.01, 0.5)
    
            l = np.array([l1, l2, l3, l4, l5])
    
            # --- final check
            C1, C2 = compute_S_10(l,k)
    
            if C1 < 0 and C2 < 0:
                return l
    
        return None    
        
    def sample_near_boundary_11(k, max_trials=1000):
    #ℓ6 + ℓ5 − 𝑘ℓ4 − (2𝑘 + 1)ℓ3 − (𝑘 + 1)ℓ2 − (𝑘 + 2)ℓ1 > 0 and ℓ1 + ℓ2 + 2ℓ3 + ℓ4 − ℓ6 > 0
        for _ in range(max_trials):
    
            # --- Step 1: sample l1,l2,l3,l4
            l1234 = np.exp(rng.standard_normal(4))
            l1, l2, l3, l4 = l1234
    
            # enforce C2 > 0
            if l1 + l2 + 2*l3 + l4 <= 1:
                continue
    
            # --- Step 2: compute lower bound for l5
            A = k*l4 + (2*k+1)*l3 + (k+1)*l2 + (k+2)*l1 - 1
    
            # choose l5
            if A <= 0:
                l5 = np.exp(rng.standard_normal()) * 0.2
            else:
                l5 = A + rng.uniform(0.01, 0.5)
    
            l = np.array([l1, l2, l3, l4, l5])
    
            # --- final check
            C1, C2 = compute_S_11(l,k)
    
            if C1 < 0 and C2 < 0:
                return l
    
        return None    
    
    def sample_near_boundary_12(k, max_trials=1000):
    
        for _ in range(max_trials):
    
            # --- Step 1: sample l1,l2,l3,l4
            l1234 = np.exp(rng.standard_normal(4))
            l1, l2, l3, l4 = l1234
    
            # enforce C2 > 0
            C2 = l1 + 2*l2 + 3*l3 + l4 - 1
            if C2 <= 0:
                continue
    
            # --- Step 2: compute lower bound for l5
            A = k*l4 + (3*k+1)*l3 + (2*k+1)*l2 + (k+1)*l1 - 1
    
            if A <= 0:
                # easy case
                l5 = np.exp(rng.standard_normal()) * 0.2
            else:
                # enforce constraint
                l5 = A + rng.uniform(0.01, 0.5)
    
            l = np.array([l1, l2, l3, l4, l5])
    
            # --- final check
            C1, C2 = compute_S_12(l,k)
    
            if C1 < 0 and C2 < 0:
                return l
    
        return None    
    
    def sample_near_boundary_13(k, max_trials=1000):
    
        for _ in range(max_trials):
    
            # --- Step 1: sample l1,l2,l3
            l1, l2, l3 = np.exp(rng.standard_normal(3))
    
            # enforce C2 > 0
            C2_base = l1 + l2 + l3
            if C2_base <= 1:
                continue
    
            # --- Step 2: choose l5 small (since subtracts in C2)
            l5 = np.exp(rng.standard_normal()) * 0.2
    
            C2 = C2_base - l5 - 1
            if C2 <= 0:
                continue
    
            # --- Step 3: solve for l4 using C1
            # C1 = 1 + 2l5 + l4 - k l3 - (k+1)l2 - (k+2)l1 > 0
    
            A = 1 + 2*l5 - k*l3 - (k+1)*l2 - (k+2)*l1
    
            # need: A + l4 > 0 → l4 > -A
            if A >= 0:
                l4 = np.exp(rng.standard_normal()) * 0.1
            else:
                l4 = -A + rng.uniform(0.01, 0.5)
    
            # assemble
            l = np.array([l1, l2, l3, l4, l5])
    
            # --- final check
            C1, C2 = compute_S_13(l,k)
    
            if C1 < 0 and C2 < 0:
                return l
    
        return None    
    
    def sample_near_boundary_14(k, max_trials=1000):
    
        for _ in range(max_trials):
    
            # --- Step 1: sample ℓ1,ℓ2,ℓ3
            l1, l2, l3 = np.exp(rng.standard_normal(3))
    
            # enforce C2 > 0
            C2_base = l1 + 2*l2 + l3
            if C2_base <= 1:
                continue
    
            # --- Step 2: choose l5 small (since it subtracts in C2)
            l5 = np.exp(rng.standard_normal()) * 0.2
    
            C2 = C2_base - l5 - 1
            if C2 <= 0:
                continue
    
            # --- Step 3: solve C1 for l4
            # C1 = 1 + 2l5 + l4 - k l3 - (2k+1)l2 - (k+1)l1 > 0
    
            A = 1 + 2*l5 - k*l3 - (2*k+1)*l2 - (k+1)*l1
    
            # need: A + l4 > 0 → l4 > -A
    
            if A >= 0:
                l4 = np.exp(rng.standard_normal()) * 0.1
            else:
                l4 = -A + rng.uniform(0.01, 0.5)
    
            l = np.array([l1, l2, l3, l4, l5])
    
            # final check
            C1, C2 = compute_S_14(l,k)
    
            if C1 < 0 and C2 < 0:
                return l
    
        return None
    
    def sample_near_boundary_15(k, max_trials=1000):
    
        for _ in range(max_trials):
            # --- Step 1: sample main contributors (control C2)
            l123 = np.exp(rng.standard_normal(3))
            l1, l2, l3 = l123
            # enforce C2 > 0 directly
            C2_base = l1 + 2*l2 + l3
            if C2_base <= 2:
                continue
            # --- Step 2: sample l5 small (since it subtracts in C2)
            l5 = np.exp(rng.standard_normal()) * 0.1
            # recompute C2
            C2 = C2_base - l5 - 2
            if C2 <= 0:
                continue
            # --- Step 3: choose l4 to satisfy C1 ≈ 0
            A = (k+1)*l1 + (2*k+1)*l2 + k*l3 - 1
            # want: C1 = A - l4 - l5 < 0 → take l4 slightly larger
            if A - l5 <= 0:
                l4 = np.exp(rng.standard_normal()) * 0.1
            else:
                l4 = A - l5 + rng.uniform(0.01, 0.5)
            l = np.array([l1, l2, l3, l4, l5])
            C1, C2 = compute_S_15(l, k)
            if C1 < 0 and C2 < 0:
                return l
        return None
    
    if(int(case)<10):
        sample_near_boundary=sample_near_boundary_1to9
    else:
        sample_near_boundary=None #depends on each case with two constraints
    
    #%%
    
          
    for k in range(kmin, kmax):    
        #print("k=",k)
        
        global_min_functional = np.inf
        global_geom = None
        best_x = None
        best_y = None
        best_ALPHA=None
        x_new = None
        y_new = None
        
        #initialization
        if(analytical==False):
            if(int(case) < 10):
                while True:
                    y0 = sample_near_boundary(k)                
                    if y0 is not None:
                        x0 = transform(y0, k)
                        if x0 is not None:
                            break
            else:
                while True:
                    if(case=="10"):
                        x0 = sample_near_boundary_10(k)
                    elif(case=="11"):
                        x0 = sample_near_boundary_11(k)
                    elif(case=="12"):
                        x0 = sample_near_boundary_12(k)
                    elif(case=="13"):
                        x0 = sample_near_boundary_13(k)
                    elif(case=="14"):
                        x0 = sample_near_boundary_14(k)
                    elif(case=="15"):
                        x0 = sample_near_boundary_15(k)                            
                    if x0 is None:
                        continue
                    # print(np.min(x0),np.max(x0))
                    x0 = np.clip(x0, xmin, xmax)
                    if(case=="10"):
                        C1, C2 = compute_S_10(x0, k)
                    elif(case=="11"):
                        C1, C2 = compute_S_11(x0, k)
                    elif(case=="12"):
                        C1, C2 = compute_S_12(x0, k)
                    elif(case=="13"):
                        C1, C2 = compute_S_13(x0, k)
                    elif(case=="14"):
                        C1, C2 = compute_S_14(x0, k)
                    elif(case=="15"):
                        C1, C2 = compute_S_15(x0, k)
                        
                    if C1 < 0 and C2 < 0:
                        break
        else:
            x0 = analytical_l_values
        
        for mc_iter in range(mc_iter_max): 
                      
            if int(case) < 10:
                if(analytical==False):
                    y_new = y0 * np.exp(step_size * rng.standard_normal(5))
                    #y_new = y0 + step_size * y0 * rng.standard_normal(5)
                    # keep inside simplex region (test even with exp)
                    if np.any(y_new < GEOM_EPS):
                        continue
                    x_new = transform(y_new, k)
                else:
                    x_new = analytical_l_values

                if(verbose>1):
                    print("range=",np.min(x_new), np.max(x_new))
                    print("S =", compute_S(x_new,k))
                if x_new is None:
                    continue
                
            else:
                if(analytical==False):
                    x_new = x0 * np.exp(step_size*rng.standard_normal(5))
                    #x_new = x0 + step_size * x0 * rng.standard_normal(5)
                else:
                    x_new = analytical_l_values

                if(verbose>1):
                    print("range=",np.min(x_new), np.max(x_new))
                if np.any(x_new < GEOM_EPS):
                    continue
                
                if(case=="10"):
                    C1, C2 = compute_S_10(x_new, k)
                elif(case=="11"):
                    C1, C2 = compute_S_11(x_new, k)
                elif(case=="12"):
                    C1, C2 = compute_S_12(x_new, k)
                elif(case=="13"):
                    C1, C2 = compute_S_13(x_new, k)
                elif(case=="14"):
                    C1, C2 = compute_S_14(x_new, k)
                elif(case=="15"):
                    C1, C2 = compute_S_15(x_new, k)
                if not (C1 < 0 and C2 < 0):
                    continue                                
        
            semiglobal_min_functional, semiglobal_geom, ALPHA = func_L(x_new, k)
        
            if semiglobal_geom is None or not np.isfinite(semiglobal_min_functional):
                continue
                   
            # --- accept if better (greedy descent) ---
            if semiglobal_min_functional < global_min_functional or global_min_functional==np.inf:
                if(analytical==False and verbose>1):
                    print("k/MC/min(L_A)/ALPHA:",k,mc_iter,step_size,semiglobal_min_functional,global_min_functional,ALPHA)                
                global_min_functional = semiglobal_min_functional
                global_geom = semiglobal_geom
                best_x = x_new.copy()
                best_y = None if y_new is None else y_new.copy()
                best_ALPHA=ALPHA.copy()
                
            if(best_x is not None):
                x0=best_x
            else:
                x0=x_new
            if(best_y is not None):
                y0=best_y
            else:
                y0=y_new
            
            if(analytical==False and global_min_functional<STOP_TOL):
                print("mc iter instability found:",mc_iter, k, global_min_functional)
                break

            if(analytical and abs(global_min_functional-analytical_Min_L)<INT_EPS):
                break    
                                                
    #%%                    
        #print("--------------------------------------------------")
        #print('case:',case)                        
        def print_constraints(x, k, case):        
            l1, l2, l3, l4, l5 = x
            S=-1
            C1=-1
            C2=-1
            if case == "1":
                S = compute_S_1(x, k)
                #print("Constraint (case 1): S =", S, "< 0")
            elif case == "2":
                S = compute_S_2(x, k)
                #print("Constraint (case 2): S =", S, "< 0")
            elif case == "3":
                S = compute_S_3(x, k)
                #print("Constraint (case 3): S =", S, "< 0")
            elif case == "4":
                S = compute_S_4(x, k)
                #print("Constraint (case 4): S =", S, "< 0")
            elif case == "5":
                S = compute_S_5(x, k)
                #print("Constraint (case 5): S =", S, "< 0")
            elif case == "6":
                S = compute_S_6(x, k)
                #print("Constraint (case 6): S =", S, "< 0")
            elif case == "7":
                S = compute_S_7(x, k)
                #print("Constraint (case 7): S =", S, "< 0")
            elif case == "8":
                S = compute_S_8(x, k)
                #print("Constraint (case 8): S =", S, "< 0")
            elif case == "9":
                S = compute_S_9(x, k)
                #print("Constraint (case 9): S =", S, "< 0")
            elif case == "10":
                C1, C2 = compute_S_10(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            elif case == "11":
                C1, C2 = compute_S_11(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            elif case == "12":
                C1, C2 = compute_S_12(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            elif case == "13":
                C1, C2 = compute_S_13(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            elif case == "14":
                C1, C2 = compute_S_14(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            elif case == "15":
                C1, C2 = compute_S_15(x, k)
                #print("Constraint C1 < 0 :", C1)
                #print("Constraint C2 < 0 :", C2)
            else:
                print("Unknown case")
    
            constraints=max(max(S,C1),C2)>0
            return constraints
            
                    
        if global_geom is not None:
    
            l1, l2, l3, l4, l5 = best_x
            if(verbose>1):                       
                print("k=",k,"l6=",1)
                print("l1=",l1)
                print("l2=",l2)
                print("l3=",l3)
                print("l4=",l4)
                print("l5=",l5)
                print("Min L=",global_min_functional)
            constraints=print_constraints(best_x, k, case)
            if(verbose>1):
                print("Degenerated:", global_geom[4])
                print("e1=", global_geom[0])
                print("e2=", global_geom[1])
                print("p1=", global_geom[2])
                print("p2=", global_geom[3])
                print("Intersect:", global_geom[6])
                print("ALPHA=",best_ALPHA)
                      
            record = {
                "case": case,
                "k": k,
                "l": {
                    "l1": format(l1, nflg),  
                    "l2": format(l2, nflg),
                    "l3": format(l3, nflg),
                    "l4": format(l4, nflg),
                    "l5": format(l5, nflg)
                },
                "min_L": format(global_min_functional, nflg),
                "constraints_ok": bool(not constraints),
                "alpha": [format(float(a), nflg) for a in best_ALPHA ], 
                "geometry": {
                    "degenerated": bool(global_geom[4]),
                    "e1": [fmt_point(P) for P in global_geom[0]],
                    "e2": [fmt_point(P) for P in global_geom[1]],
                    "p1": fmt_point(global_geom[2]),
                    "p2": fmt_point(global_geom[3]),
                    "intersections": [fmt_point(P) for P in global_geom[6]]                    
                }
            }
            
            if(analytical==False):
                filename = "Delzant_solution_case_"+str(case)+".json"
                
                if k>kmin and os.path.exists(filename):
                    with open(output_dir / filename, "r") as f:
                        data = json.load(f)
                
                    # Convert old single-record format to list format
                    if isinstance(data, dict):
                        data = [data]
                else:
                    data = []
                
                # Add new solution
                data.append(record)
                
                with open(output_dir / filename, "w") as f:
                    json.dump(data, f, indent=4)
            
        #print("--------------------------------------------------")
        
        if(global_min_functional<STOP_TOL):
            break
    
    
    #%%
    #synthesis for all k
    if(False and analytical==False):
        filename = "Delzant_solution_case_"+str(case)+".json"
        with open(output_dir / filename, "r") as f:                                      
            data = json.load(f)
        
        df = pd.DataFrame(data)
        
        #print(df[["case", "k", "min_L", "constraints_ok"]])
        
        fig, ax = plt.subplots(figsize=(10, 5))
        
        # Main curve
        ax.plot(df["k"], df["min_L"].astype(float), "o-", lw=2, label="min_L")
        
        # Highlight invalid solutions
        bad = ~df["constraints_ok"]
        ax.scatter(
            df.loc[bad, "k"],
            df.loc[bad, "min_L"],
            color="red",
            s=80,
            label="constraints failed",
            zorder=3,
        )
        
        # Highlight optimum
        best_idx = df["min_L"].idxmin()
        ax.scatter(
            df.loc[best_idx, "k"],
            float(df.loc[best_idx, "min_L"]),
            color="gold",
            edgecolor="black",
            s=150,
            label="best solution",
            zorder=4,
        )
        
        #ax.set_yscale("log")  # logarithmic scale (base 10)        
        df["min_L"] = pd.to_numeric(df["min_L"])
        best = df["min_L"].min()        
        ax.set_ylim(best-3*abs(best), 10*abs(float(best))  )
        ax.set_xlabel("k")
        ax.set_ylabel("min_L")
        ax.set_title(f"Case {case}: objective value vs k")
        ax.grid(True, alpha=0.3)
        ax.legend()
        
        plt.tight_layout()
        plt.show()
    
    #%%    
    #saving most unstable for a posteriori test with analytical==True
       
    if(analytical==False):
        # Best solution
        #best = min(data, key=lambda x: float(x["min_L"]))    
        best = min(
        (
            r for r in data
            if (not r["geometry"]["degenerated"])
            and r["constraints_ok"]
        ),
        key=lambda x: float(x["min_L"])
        )       
    
        if(verbose>1):
            print("Best unstable solution:")    
            print(f"k = {best['k']}")
            print(f"min_L = {best['min_L']}")
            print(f"constraints_ok = {best['constraints_ok']}")
            
            print("\nl:")
            for name, value in best["l"].items():
                print(f"  {name}: {value}")
            
            print("\ngeometry:")
            print(f"  degenerated: {best['geometry']['degenerated']}")
            print(f"  e1: {best['geometry']['e1']}")
            print(f"  e2: {best['geometry']['e2']}")
            print(f"  p1: {best['geometry']['p1']}")
            print(f"  p2: {best['geometry']['p2']}")
            print(f"  intersections: {best['geometry']['intersections']}")
            
        best_k=best['k']
            
        filename = "Delzant_solution_case_"+str(case)+"_k"+str(best_k)+"_best.json"    
        with open(output_dir / filename, "w") as f:
            json.dump(best, f, indent=4)
        
    else:
        #read _best file 
        filename=f"Delzant_solution_case_{case}_k"+str(analytical_k_values)+"_best.json"
        print("read file for a posteriori testing:",filename)       
        with open(output_dir / filename, "r") as f:
            best = json.load(f)
    
    #%%
    #a posteriori test with sympy
    
    # =====================================================
    # TEST BEST
    # =====================================================
    
    data = best
        
    k = Integer(data["k"])
    
    min_L = Float(data["min_L"], nfl)
    
    l1 = Float(data["l"]["l1"], nfl)
    l2 = Float(data["l"]["l2"], nfl)
    l3 = Float(data["l"]["l3"], nfl)
    l4 = Float(data["l"]["l4"], nfl)
    l5 = Float(data["l"]["l5"], nfl)
    
    l6 = Float(1, nfl)
    
    alpha_json = [Float(a, nfl) for a in data["alpha"]]
    
    # =====================================================
    # REBUILD POLYGON
    # =====================================================
    
    x, y = symbols("x y")
    
    verts_np = build_vertices(
        float(l1),
        float(l2),
        float(l3),
        float(l4),
        float(l5),
        float(l6),
        int(k)
    )
    
    verts = [
        (Float(v[0], nfl), Float(v[1], nfl))
        for v in verts_np
    ]        
    
    x0,y0 = verts[0]
    x1,y1 = verts[1]
    x2,y2 = verts[2]
    x3,y3 = verts[3]
    x4,y4 = verts[4]
    x5,y5 = verts[5]
    x6,y6 = verts[6]
    x7,y7 = verts[7]
    
    # =====================================================
    # NUMERICAL CHECK
    # =====================================================
    
    alpha_np = np.array(
        [float(a) for a in data["alpha"]],
        dtype=float
    )
    
    p1_np = np.array(data["geometry"]["p1"], dtype=float)
    p2_np = np.array(data["geometry"]["p2"], dtype=float)
    
    LA_num, B_num, A_num, _ = compute_LA_global(
        p1_np,
        p2_np,
        alpha=alpha_np
    )
    
    if(verbose>1):
        print("\n--- NUMERICAL CHECK ---")
        print("Numeric min_L =", data["min_L"])
        print("Recomputed = {:.60e}".format(LA_num))
        print("Difference =", abs(Float(str(LA_num), nfl) - Float(data["min_L"], nfl)
    )
    )
    
    # =====================================================
    # STORED LINE
    # =====================================================
    
    p1 = data["geometry"]["p1"]
    p2 = data["geometry"]["p2"]
    
    px1 = Float(p1[0], nfl)
    py1 = Float(p1[1], nfl)
    
    px2 = Float(p2[0], nfl)
    py2 = Float(p2[1], nfl)
    
    aL = py2 - py1
    bL = -(px2 - px1)
    cL = px2*py1 - py2*px1
    
    Lexpr = expand(aL*x + bL*y + cL)
    
    Aexpr = expand(
        alpha_json[0]
        + alpha_json[1]*x
        + alpha_json[2]*y
    )
    
    # =====================================================
    # BUILD P+
    # =====================================================
    
    def build_positive_polygon():
    
        Pplus = []
    
        n = len(verts)
    
        for i in range(n):
    
            A = verts[i]
            B = verts[(i+1)%n]
    
            sA = float(
                N(
                    Lexpr.subs({
                        x:A[0],
                        y:A[1]
                    }),
                    nfl
                )
            )
    
            sB = float(
                N(
                    Lexpr.subs({
                        x:B[0],
                        y:B[1]
                    }),
                    nfl
                )
            )
    
            if sA >= 0:
                Pplus.append(A)
    
            if sA*sB < 0:
    
                u = symbols(f"u_{i}")
    
                dx = B[0]-A[0]
                dy = B[1]-A[1]
    
                root = solve(
                    expand(
                        Lexpr.subs({
                            x:A[0] + u*dx,
                            y:A[1] + u*dy
                        })
                    ),
                    u
                )[0]
    
                I = (
                    simplify(A[0] + root*dx),
                    simplify(A[1] + root*dy)
                )
    
                Pplus.append(I)
    
        return Pplus
    
    Pplus = build_positive_polygon()
    
    #print("\nPplus vertices =", len(Pplus))
    
    # =====================================================
    # TRIANGLE INTEGRAL
    # =====================================================
    
    def triangle_integral(expr,P,Q,R):
    
        s,t = symbols("s t")
    
        X = P[0] + s*(Q[0]-P[0]) + t*(R[0]-P[0])
        Y = P[1] + s*(Q[1]-P[1]) + t*(R[1]-P[1])
    
        J = Abs(
            Matrix([
                [Q[0]-P[0], R[0]-P[0]],
                [Q[1]-P[1], R[1]-P[1]]
            ]).det()
        )
    
        expr2 = expand(expr.subs({x:X,y:Y}))
    
        return integrate(
            integrate(
                expr2*J,
                (t,0,1-s)
            ),
            (s,0,1)
        )
    
    # =====================================================
    # EXACT SYMBOLIC L_A
    # =====================================================
    
    def exact_LA(Pplus):
    
        integrand = expand(Aexpr*Lexpr)
    
        area_term = 0
    
        P0 = Pplus[0]
    
        for i in range(1,len(Pplus)-1):
    
            area_term += triangle_integral(
                integrand,
                P0,
                Pplus[i],
                Pplus[i+1]
            )
    
        primitive_dirs = get_primitive_dirs(case, k)
    
        verts_closed = verts + [verts[0]]
    
        boundary_term = 0
    
        for e in range(8):
    
            Aedge = verts_closed[e]
            Bedge = verts_closed[e+1]
    
            vx,vy = primitive_dirs[e]
    
            dx = Bedge[0]-Aedge[0]
            dy = Bedge[1]-Aedge[1]
    
            sigma = abs(dx/vx) if abs(vx)>0 else abs(dy/vy)
    
            t = symbols(f"tt_{e}")
    
            X = Aedge[0] + t*dx
            Y = Aedge[1] + t*dy
    
            Ledge = expand(
                Lexpr.subs({x:X,y:Y})
            )
    
            sA = float(N(Lexpr.subs({x:Aedge[0],y:Aedge[1]}),nfl))
            sB = float(N(Lexpr.subs({x:Bedge[0],y:Bedge[1]}),nfl))
    
            roots = solve(Ledge,t)
    
            if sA > 0 and sB > 0:
    
                boundary_term += sigma*integrate(
                    Ledge,
                    (t,0,1)
                )
    
            elif len(roots)==1:
    
                r = roots[0]
    
                rv = float(N(r))
                if 0.0 < rv < 1.0:
    
                    if sA > 0:
    
                        boundary_term += sigma*integrate(
                            Ledge,
                            (t,0,r)
                        )
    
                    elif sB > 0:
    
                        boundary_term += sigma*integrate(
                            Ledge,
                            (t,r,1)
                        )
    
        return (
            simplify(boundary_term-area_term),
            simplify(boundary_term),
            simplify(area_term)
        )
    
    LA_sym, B_sym, A_sym = exact_LA(Pplus)
    
    LA_numeric = Float(data["min_L"], nfl)
    
    # print("\n--- SYMBOLIC CHECK ---")
    # print("Boundary =", N(B_sym,nfl))
    # print("Area     =", N(A_sym,nfl))
    print("Symbolic L_A      =", float(N(LA_sym,nfl)))    
    print("Numeric L_A     =", float(LA_numeric))
    rel_err = abs(LA_sym - LA_numeric) / (1.e-30 + abs(LA_numeric))
    print("Relative error =", float(N(rel_err,nfl)))
    best_k=best['k']
    best_constraints_ok=best['constraints_ok']
        
    return LA_numeric, rel_err,  best_k, best_constraints_ok


#%%

if __name__ == "__main__":
    """
    either parameteric search for k for given case, or try all cases at given k=kmin, set kmax=kmin+1:
    for multi-case trials freeze kmin and kmax=kmin+1
    """
    
    cases=np.arange(1,16)  #cases during optimization or testing with analytical=True

    kmin=20
    kmax=kmin+1
    
    analytical=False
    #if analytical: 
    #kmin will be set to analytical_k_values found previously read in load_analytical_solution
    #and kmax will be set to analytical_k_values+1

    res_cases=[]
          
    for icase in cases: 
        case=str(icase)
        print("case:",case)
        
        LA_numeric, rel_err,  best_k, constraints_ok = octogone_donaldson_delzant_fct(analytical, case, kmin, kmax)
        res_cases.append([icase, best_k, float(LA_numeric), float(rel_err), constraints_ok])

    res_cases=np.array(res_cases)        
    y1 = res_cases[:,1]
    y2 = res_cases[:,2]
    
    fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharex=True)
    
    ax[0].plot(res_cases[:,0], y1, "o-", label="k min_LA")
    ax[0].set_xlabel("cases")
    ax[0].legend()
    
    ax[1].plot(res_cases[:,0], y2, "o-", label="min_LA")
    ax[1].set_xlabel("cases")
    ax[1].legend()
    
    plt.tight_layout()
    plt.show()