import numpy as np
import matplotlib.pyplot as plt
import sympy

def generate_primes(k):
    primes = []
    num = 2
    while len(primes) < k:
        if all(num % p != 0 for p in primes if p*p <= num):
            primes.append(num)
        num += 1
    return primes

def wave_function(x, primes, amplitude_scale=0.25, invert=True):
    total = np.zeros_like(x, dtype=float)
    for p in primes:
        min_amp = p / 1000
        amplitude = min_amp + amplitude_scale * (p - min_amp)
        wave = amplitude * np.sin(2 * np.pi * x / p)
        if invert:
            wave = np.where(wave > 0, -wave, wave)
        total += wave
    return total

def find_pattern_primes(k):
    primes = generate_primes(k)
    N_k = np.prod(primes)
    pattern_primes = []
    p = primes[-1] + 1
    while p <= N_k:
        if all(p % prime != 0 for prime in primes):
            pattern_primes.append(p)
        p = sympy.nextprime(p)
    return primes, N_k, pattern_primes

def check_symmetry(k):
    primes = generate_primes(k)
    N_k = np.prod(primes)
    center = N_k / 2
    _, _, pattern_primes = find_pattern_primes(k)
    symmetry_holds = True
    violations = []
    for p_val in pattern_primes:
        if p_val > center:
            continue
        mirror = N_k - p_val
        if mirror <= N_k and mirror not in pattern_primes:
            symmetry_holds = False
            violations.append((p_val, mirror))
    return symmetry_holds, violations, pattern_primes

def plot_wave_function(k, primes, N_k):
    num_points = 5000
    x = np.linspace(0.4 * N_k, 0.6 * N_k, num_points)
    y = wave_function(x, primes)
    
    plt.figure(figsize=(14, 8))
    plt.plot(x, y, 'b-', linewidth=1.5, label=f'F_k(x) для k={k}')
    plt.axvline(N_k/2, color='r', linestyle='--', label='Центр симметрии')
    
    # Отметить простые числа
    _, _, pattern_primes = find_pattern_primes(k)
    for p in pattern_primes:
        if 0.4 * N_k <= p <= 0.6 * N_k:
            plt.axvline(p, color='g', alpha=0.4, linestyle=':')
    
    plt.title(f'Волновая функция для первых {k} простых чисел (N_k={N_k})', fontsize=14)
    plt.xlabel('x', fontsize=12)
    plt.ylabel('F_k(x)', fontsize=12)
    plt.legend()
    plt.grid(True, linestyle='--', alpha=0.5)
    plt.savefig(f'wave_function_k={k}.png', dpi=300)
    plt.show()

if __name__ == "__main__":
    # Проверка симметрии для различных k
    for k in [3, 4, 5, 6]:
        symmetry_holds, violations, pattern_primes = check_symmetry(k)
        primes = generate_primes(k)
        N_k = np.prod(primes)
        
        print(f"\nАнализ для k={k} (N_k={N_k})")
        print(f"Количество простых в паттерне: {len(pattern_primes)}")
        
        if symmetry_holds:
            print("Симметрия подтверждена: все числа имеют зеркальные пары")
        else:
            print(f"Нарушения симметрии: {len(violations)}")
            print("Примеры нарушений:")
            for i, (p_val, mirror) in enumerate(violations[:3]):
                print(f"  {p_val} → {mirror} (простое: {sympy.isprime(mirror)})")
        
        # Визуализация волновой функции
        plot_wave_function(k, primes, N_k)
    
    # Анализ сходимости
    convergence_results = []
    k_values = [3, 4, 5, 6, 7, 8, 9, 10]
    
    for k in k_values:
        primes = generate_primes(k)
        N_k = np.prod(primes)
        convergence_results.append(N_k)
        print(f"k={k}: N_k={N_k}")
    
    # Визуализация роста N_k
    plt.figure(figsize=(10, 6))
    plt.plot(k_values, convergence_results, 'bo-', markersize=8)
    plt.yscale('log')
    plt.title('Экспоненциальный рост N_k', fontsize=14)
    plt.xlabel('Количество простых чисел (k)', fontsize=12)
    plt.ylabel('N_k (log scale)', fontsize=12)
    plt.grid(True, which="both", linestyle='--', alpha=0.5)
    
    for k, N_k in zip(k_values, convergence_results):
        plt.annotate(f"{N_k:.1e}", (k, N_k), textcoords="offset points", xytext=(0,10), ha='center')
    
    plt.tight_layout()
    plt.savefig('nk_growth.png', dpi=300)
    plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress
from math import log

def generate_primes(k):
    primes = []
    num = 2
    while len(primes) < k:
        if all(num % p != 0 for p in primes if p*p <= num):
            primes.append(num)
        num += 1
    return primes

def wave_function(x, primes, amplitude_scale=0.25, invert=True):
    total = np.zeros_like(x, dtype=float)
    for p in primes:
        min_amp = p / 1000
        amplitude = min_amp + amplitude_scale * (p - min_amp)
        wave = amplitude * np.sin(2 * np.pi * x / p)
        if invert:
            wave = np.where(wave > 0, -wave, wave)
        total += wave
    return total

def fractal_dimension_box_counting(x, y, num_scales=12):
    x_norm = (x - np.min(x)) / (np.max(x) - np.min(x))
    y_norm = (y - np.min(y)) / (np.max(y) - np.min(y))
    points = np.column_stack((x_norm, y_norm))

    scales = np.logspace(-2, 0, num_scales, endpoint=False)[::-1]
    counts = []
    for scale in scales:
        grid = set()
        for point in points:
            cell_x = int(point[0] / scale)
            cell_y = int(point[1] / scale)
            grid.add((cell_x, cell_y))
        counts.append(len(grid))

    log_scales = np.log(1 / scales)
    log_counts = np.log(counts)

    slope, intercept, r_value, _, _ = linregress(log_scales, log_counts)
    return slope, r_value**2, log_scales, log_counts

if __name__ == "__main__":
    k_values = [3, 4, 5, 6, 7, 8]
    results = []

    # Создаем фигуру с двумя подграфиками
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), gridspec_kw={'height_ratios': [2, 1]})

    for k in k_values:
        primes = generate_primes(k)
        N_k = np.prod(primes)

        num_points = min(100000, max(20000, int(50000 / log(k+1))))
        start_frac = 0.45 - 0.01 * k
        end_frac = 0.55 + 0.01 * k

        x = np.linspace(start_frac * N_k, end_frac * N_k, num_points)
        y = wave_function(x, primes)

        dim, r2, log_scales, log_counts = fractal_dimension_box_counting(x, y)

        results.append({'k': k, 'N_k': N_k, 'dim': dim, 'r2': r2})

        ax1.plot(k, dim, 'o', label=f'k={k}' if k in [3,5,7,8] else None)

        if k == 7:
            inset_ax = ax1.inset_axes([0.6, 0.15, 0.35, 0.35])
            inset_ax.plot(log_scales, log_counts, 'bo-', label='Данные')
            fit_line = dim * log_scales + (log_counts[0] - dim * log_scales[0])
            inset_ax.plot(log_scales, fit_line, 'r--', label=f'Линейная регрессия\nD={dim:.4f}\nR²={r2:.4f}')
            inset_ax.set_xlabel('log(1/scale)', fontsize=8)
            inset_ax.set_ylabel('log(N)', fontsize=8)
            inset_ax.tick_params(axis='both', which='major', labelsize=8)
            inset_ax.grid(True, linestyle='--', linewidth=0.5)
            inset_ax.legend(fontsize=7)
            inset_ax.set_title('Логарифмическая зависимость (k=7)', fontsize=9)

        print(f"k={k}: N_k={N_k}, D={dim:.4f}, R²={r2:.4f}")

    ax1.axhline(y=1.73, color='r', linestyle='--', label='Целевое значение 1.73')
    ax1.set_title('Фрактальная размерность волновых паттернов', fontsize=16)
    ax1.set_xlabel('Количество простых чисел (k)', fontsize=14)
    ax1.set_ylabel('Фрактальная размерность (D)', fontsize=14)
    ax1.grid(True, linestyle='--', alpha=0.5)
    ax1.legend(fontsize=10)
    ax1.set_xticks(k_values)
    ax1.set_ylim(0.9, 2.1)

    # Волновая функция для k=7
    k7 = next(res for res in results if res['k'] == 7)
    primes = generate_primes(7)
    N_k = k7['N_k']
    x7 = np.linspace(0.45 * N_k, 0.55 * N_k, 50000)
    y7 = wave_function(x7, primes)

    ax2.plot(x7, y7, color='blue', linewidth=0.8)
    ax2.set_title(f'Волновая функция для k=7 (N_k={N_k})', fontsize=16)
    ax2.set_xlabel('Позиция', fontsize=14)
    ax2.set_ylabel('Амплитуда', fontsize=14)
    ax2.grid(True, linestyle='--', alpha=0.5)

    textstr = f'Фрактальная размерность: {k7["dim"]:.4f}\nR²: {k7["r2"]:.4f}'
    props = dict(boxstyle='round', facecolor='white', alpha=0.85, edgecolor='gray')
    ax2.text(0.5, 0.05, textstr, transform=ax2.transAxes, fontsize=12,
             verticalalignment='bottom', horizontalalignment='center', bbox=props)

    plt.tight_layout()
    plt.show()

    # Средняя размерность для k≥5
    large_k_dims = [res['dim'] for res in results if res['k'] >= 5]
    mean_dim = np.mean(large_k_dims)
    print(f"\nСредняя фрактальная размерность для k≥5: {mean_dim:.4f}")
    print(f"Отклонение от 1.73: {abs(mean_dim - 1.73):.4f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import wasserstein_distance, ks_2samp, gaussian_kde, entropy
from scipy.optimize import curve_fit
from sklearn.utils import resample
import sympy
import time
import numba
from tqdm import tqdm
import logging

# Настройка логирования
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger()

# Ускорение вычислений
@numba.njit
def generate_primes(k):
    primes = []
    num = 2
    while len(primes) < k:
        is_prime = True
        for p in primes:
            if p * p > num:
                break
            if num % p == 0:
                is_prime = False
                break
        if is_prime:
            primes.append(num)
        num += 1
    return primes

def get_high_precision_data(k, base_sample_size=100000):
    """Улучшенный сбор данных с адаптивным размером выборки"""
    primes = generate_primes(k)
    last_prime = primes[-1]
    N_k = np.prod(primes)
    
    # Адаптивный размер выборки
    sample_size = min(10**7, base_sample_size * (k//3 + 1))
    max_size = last_prime + sample_size
    
    # Для больших k расширяем диапазон
    if k > 10:
        max_size = min(N_k, max_size * 3)
    
    logger.debug(f"k={k}: last_prime={last_prime}, sample_size={sample_size}, max_size={max_size}")
    
    # Генерация простых чисел
    pattern_primes = list(sympy.primerange(last_prime + 1, max_size))
    logger.debug(f"Найдено {len(pattern_primes)} простых чисел")
    
    # Определение количества частей
    min_points_per_section = 50
    num_sections = min(3, len(pattern_primes) // min_points_per_section)
    
    if num_sections < 2:
        logger.warning(f"Недостаточно данных для k={k}: {len(pattern_primes)} простых чисел")
        return None
    
    logger.info(f"k={k}: {len(pattern_primes)} простых, разделение на {num_sections} части")
    
    # Разделение данных
    section_size = len(pattern_primes) // num_sections
    intervals = []
    
    for i in range(num_sections):
        start_idx = i * section_size
        end_idx = (i+1) * section_size if i < num_sections-1 else len(pattern_primes)
        section = pattern_primes[start_idx:end_idx]
        intervals.append(np.diff(section))
    
    return intervals

def js_divergence(p, q):
    """Расчет дивергенции Дженсена-Шеннона"""
    p = np.asarray(p)
    q = np.asarray(q)
    p = p / p.sum()
    q = q / q.sum()
    m = 0.5 * (p + q)
    return 0.5 * (entropy(p, m) + entropy(q, m))

def calculate_advanced_metrics(intervals):
    """Расчет метрик с защитой от ошибок"""
    results = {}
    num_sections = len(intervals)
    
    # Создаем пары для сравнения
    pairs = []
    for i in range(num_sections):
        for j in range(i+1, num_sections):
            pairs.append((i, j))
    
    # Объединение всех данных
    all_intervals = np.concatenate(intervals)
    hist_range = (0, np.max(all_intervals)*1.2)
    bins = np.linspace(hist_range[0], hist_range[1], 100)
    
    for i, j in pairs:
        key = f"comp_{i}_{j}"
        try:
            dist1, dist2 = intervals[i], intervals[j]
            
            # Пропуск пустых распределений
            if len(dist1) < 10 or len(dist2) < 10:
                continue
                
            # Нормализация
            norm1 = dist1 / np.mean(dist1)
            norm2 = dist2 / np.mean(dist2)
            
            # Расстояния
            results[f"{key}_emd"] = wasserstein_distance(norm1, norm2)
            results[f"{key}_ks"] = ks_2samp(norm1, norm2).statistic
            
            # JS-дивергенция
            hist1, _ = np.histogram(dist1, bins=bins, range=hist_range)
            hist2, _ = np.histogram(dist2, bins=bins, range=hist_range)
            results[f"{key}_js"] = js_divergence(hist1, hist2)
        except Exception as e:
            logger.warning(f"Ошибка расчета для {key}: {str(e)}")
    
    # Коэффициенты вариации
    for i in range(num_sections):
        try:
            results[f"cv_{i}"] = np.std(intervals[i]) / np.mean(intervals[i])
        except:
            results[f"cv_{i}"] = np.nan
            
    results["overall_cv"] = np.std(all_intervals) / np.mean(all_intervals)
    
    # Оценка плотности (только если есть хотя бы 2 распределения)
    if num_sections >= 2:
        try:
            grid = np.linspace(0, np.max(all_intervals)*1.2, 500)
            pdfs = []
            for dist in intervals:
                if len(dist) > 10:
                    kde = gaussian_kde(dist, bw_method='scott')
                    pdfs.append(kde(grid))
            
            if len(pdfs) >= 2:
                pdf_diff = []
                for i in range(len(pdfs)):
                    for j in range(i+1, len(pdfs)):
                        pdf_diff.append(np.mean(np.abs(pdfs[i]-pdfs[j])))
                results["pdf_diff"] = np.mean(pdf_diff)
        except:
            results["pdf_diff"] = np.nan
    
    return results

def bootstrap_ci(data, metric_func, n_bootstrap=1000, ci=95):
    """Бутстрэп доверительных интервалов"""
    stats = []
    n = len(data)
    
    for _ in range(n_bootstrap):
        sample = resample(data, replace=True, n_samples=n)
        stats.append(metric_func(sample))
    
    alpha = (100 - ci) / 2
    lower = np.percentile(stats, alpha)
    upper = np.percentile(stats, 100 - alpha)
    return np.mean(stats), lower, upper

def asymptotic_model(k, a, b, c):
    """Экспоненциальная модель"""
    return a * np.exp(-b * k) + c

def analyze_quantum_precision(max_k=20, precision_samples=10):
    """Анализ квантового предела с улучшенной логикой"""
    results = []
    k_values = list(range(5, max_k + 1))  # Начинаем с k=5
    
    logger.info("Запуск высокоточного анализа...")
    start_time = time.time()
    
    # Для кросс-валидации
    cross_val_metrics = {'emd': [], 'js': [], 'pdf_diff': []}
    
    for k in tqdm(k_values, desc="Анализ k"):
        interval_data = get_high_precision_data(k)
        
        if interval_data is None:
            logger.warning(f"Пропуск k={k}: недостаточно данных")
            continue
            
        logger.debug(f"Получено {len(interval_data)} интервалов")
            
        # Многократное вычисление метрик
        metrics_sum = {}
        metrics_count = {}
        
        for _ in range(precision_samples):
            try:
                metrics = calculate_advanced_metrics(interval_data)
                for key, value in metrics.items():
                    if np.isfinite(value):
                        metrics_sum[key] = metrics_sum.get(key, 0) + value
                        metrics_count[key] = metrics_count.get(key, 0) + 1
            except Exception as e:
                logger.error(f"Ошибка расчета метрик для k={k}: {str(e)}")
        
        # Усреднение результатов
        avg_metrics = {}
        for key in metrics_sum:
            avg_metrics[key] = metrics_sum[key] / metrics_count[key]
        
        avg_metrics['k'] = k
        results.append(avg_metrics)
        
        # Кросс-валидация (только если есть несколько интервалов)
        if len(interval_data) >= 2:
            try:
                # Используем первые 2/3 как train, последнюю 1/3 как test
                split_idx = int(len(interval_data) * 2 / 3)
                train_data = interval_data[:split_idx]
                test_data = interval_data[split_idx:]
                
                # Вычисляем метрики
                train_metrics = calculate_advanced_metrics(train_data)
                test_metrics = calculate_advanced_metrics(test_data)
                
                # Сохраняем разницу
                for m in ['emd', 'js', 'pdf_diff']:
                    for key in train_metrics:
                        if m in key:
                            diff = abs(train_metrics.get(key, 0) - test_metrics.get(key, 0))
                            cross_val_metrics[m].append(diff)
                            break
            except Exception as e:
                logger.error(f"Ошибка кросс-валидации: {str(e)}")
    
    logger.info(f"Анализ завершен за {time.time() - start_time:.2f} сек")
    
    if not results:
        logger.error("Нет данных для анализа")
        return
    
    # Анализ стабилизации
    k_values = [r['k'] for r in results]
    
    # Основная метрика - используем первую доступную EMD
    emd_key = None
    for key in results[0]:
        if "emd" in key:
            emd_key = key
            break
    
    if emd_key:
        emd_values = [r[emd_key] for r in results]
    else:
        logger.error("EMD метрики не найдены")
        return
    
    # Улучшенный критерий стабилизации
    stability_threshold = 0.01
    stable_zone_start = None
    
    for i in range(3, len(emd_values)):
        # Проверяем стабильность в окне из 3 точек
        window = emd_values[i-2:i+1]
        if max(window) - min(window) < stability_threshold:
            stable_zone_start = k_values[i-2]
            break
    
    if stable_zone_start is None:
        logger.warning("Стабилизация не обнаружена")
        stable_zone = []
    else:
        stable_zone = [r for r in results if r['k'] >= stable_zone_start]
        logger.info(f"Зона стабилизации начинается с k={stable_zone_start}")
        
        # Расчет квантового предела
        try:
            emd_stable_vals = [r[emd_key] for r in stable_zone]
            quantum_limit, lower_ci, upper_ci = bootstrap_ci(
                emd_stable_vals,
                np.mean
            )
            
            logger.info("\nТОЧНЫЙ КВАНТОВЫЙ ПРЕДЕЛ ТОЧНОСТИ")
            logger.info(f"Величина предела: {quantum_limit:.6f}")
            logger.info(f"95% доверительный интервал: [{lower_ci:.6f}, {upper_ci:.6f}]")
            logger.info(f"Точность оценки: ±{(upper_ci - lower_ci)/2:.6f}")
            
            # Сохранение результатов
            with open("quantum_limit_results.txt", "w") as f:
                f.write(f"Quantum Precision Limit = {quantum_limit:.8f} ± {(upper_ci - lower_ci)/2:.8f}\n")
                f.write(f"Confidence Interval: [{lower_ci:.8f}, {upper_ci:.8f}]\n")
                f.write(f"Stabilization starts at k={stable_zone_start}\n")
                f.write(f"Based on {len(stable_zone)} data points\n")
        except Exception as e:
            logger.error(f"Ошибка расчета предела: {str(e)}")
    
    # Визуализация
    plt.figure(figsize=(15, 10))
    
    # Основной график EMD
    plt.subplot(2, 2, 1)
    plt.plot(k_values, emd_values, 'bo-', label='EMD')
    if stable_zone and 'quantum_limit' in locals():
        plt.axhline(y=quantum_limit, color='r', linestyle='--', label='Квантовый предел')
        plt.fill_between(k_values, lower_ci, upper_ci, color='r', alpha=0.1)
    plt.title('Эволюция расстояния EMD', fontsize=14)
    plt.xlabel('k', fontsize=12)
    plt.ylabel('EMD', fontsize=12)
    plt.grid(True)
    plt.legend()
    
    # График JS-дивергенции
    plt.subplot(2, 2, 2)
    js_values = [r.get('comp_0_1_js', np.nan) for r in results]
    plt.plot(k_values, js_values, 'co-')
    plt.title('JS-дивергенция между секциями', fontsize=14)
    plt.xlabel('k', fontsize=12)
    plt.ylabel('JS-дивергенция', fontsize=12)
    plt.grid(True)
    
    # Экспоненциальная аппроксимация
    plt.subplot(2, 2, 3)
    if len(emd_values) > 4:
        try:
            k_arr = np.array(k_values)
            emd_arr = np.array(emd_values)
            
            # Начальные параметры
            a_guess = emd_arr[0] - emd_arr[-1]
            b_guess = 0.2
            c_guess = emd_arr[-1]
            
            params, _ = curve_fit(asymptotic_model, k_arr, emd_arr, 
                                 p0=[a_guess, b_guess, c_guess], maxfev=10000)
            a, b, c = params
            quantum_limit_exp = c
            
            # Прогноз
            k_ext = np.linspace(min(k_values), max(k_values) + 5, 100)
            emd_fit = asymptotic_model(k_ext, a, b, c)
            
            plt.plot(k_values, emd_values, 'bo-', label='Данные')
            plt.plot(k_ext, emd_fit, 'r--', label=f'Модель: {a:.3f}·exp(-{b:.3f}k) + {c:.3f}')
            plt.axhline(y=c, color='g', linestyle='-.', label=f'Асимптота: {c:.5f}')
            plt.title('Экспоненциальная аппроксимация', fontsize=14)
            plt.xlabel('k', fontsize=12)
            plt.ylabel('EMD', fontsize=12)
            plt.grid(True)
            plt.legend()
        except Exception as e:
            logger.error(f"Ошибка аппроксимации: {str(e)}")
            plt.plot(k_values, emd_values, 'bo-')
            plt.title('EMD (без аппроксимации)', fontsize=14)
            plt.grid(True)
    
    plt.tight_layout()
    plt.savefig('quantum_limit_results.png', dpi=300)
    plt.show()

if __name__ == "__main__":
    # Для отладки
    logger.setLevel(logging.INFO)
    
    analyze_quantum_precision(max_k=20, precision_samples=5)