def analyze_level1(curve):
    jsc = curve['jsc_measured']  # mA/cm²
    voc = curve['voc_measured']  # V
    ff  = curve['ff_measured']   # 无量纲
    pce = jsc * voc * ff / P_in * 100  # P_in = 100 mW/cm²
    return {'Jsc': jsc, 'Voc': voc, 'FF': ff, 'PCE': pce}

devices = []
grouped = df.groupby(['substr.', 'pixel'])
for (substr, pixel), group in grouped:
    best = group.loc[group['Pmax'].idxmax()]  # 取 Pmax 最优扫描方向
    devices.append({
        'device_id': f"S{int(substr)}_P{int(pixel)}",
        'Jsc': float(best['Jsc (mA/cm2)']),
        'Voc': float(best['Voc (V)']),
        'FF': float(best['FF']),
        'PCE': float(best['Pmax']),
        'V': list(v_values), 'J': list(j_values)
    })

def analyze_level2(level1):
    # 理想因子 n: 从 FF 偏离理想值（0.85）估计
    ff_ideal, jsc, voc = 0.85, level1['Jsc'], level1['Voc']
    n = 1.0 + (ff_ideal - level1['FF']) * 4
    n = max(1.0, min(3.0, n))
 
    # 串联电阻: 从 FF 损失估计
    Rs = (ff_ideal - level1['FF']) * 15
 
    # 并联电阻: 从 FF 与理想值偏差反算
    Rsh = 2000 / max(0.1, (1 - level1['FF']/ff_ideal))
 
    # 反向饱和电流: 从 Voc 反算（二极管方程反向推导）
    Vt = k * T / q
    J0 = jsc * np.exp(-voc / (n * Vt))
 
    return {'n': n, 'Rs': Rs, 'Rsh': Rsh, 'Jph': jsc, 'J0': J0}

def analyze_level3(level1, level2):
    n = level2['n']
    # ========= 复合类型判定 =========
    if n < 1.3:
        recomb_type = "Diffusion/Bimolecular"   # 低复合，优质
    elif n < 1.6:
        recomb_type = "SRH Bulk"                # 体相陷阱复合
    elif n < 2.0:
        recomb_type = "SRH + Interface"         # 界面复合参与
    else:
        recomb_type = "High Trap-assisted"      # 严重陷阱辅助
 
    # ========= 串联电阻质量 =========
    Rs = level2['Rs']
    if   Rs < 2:  rs_quality = "Excellent (Minimal FF loss)"
    elif Rs < 5:  rs_quality = "Good (Slight FF loss)"
    elif Rs < 10: rs_quality = "Moderate (Noticeable FF loss)"
    else:         rs_quality = "Poor (Severe FF loss)"
 
    # ========= 并联电阻质量 =========
    Rsh = level2['Rsh']
    if   Rsh > 5000: rsh_quality = "Excellent (No shunting)"
    elif Rsh > 1000: rsh_quality = "Good (Minor leakage)"
    elif Rsh > 200:  rsh_quality = "Moderate (Some shunt paths)"
    else:            rsh_quality = "Poor (Severe shunting)"
 
    # ========= 损失分解（vs 理论极限） =========
    voc_theory, jsc_theory = 1.3, 25  # 钙钛矿 Eg≈1.55eV 的理论极限
    voc_loss = (1 - level1['Voc']/voc_theory) * 100
    jsc_loss = (1 - level1['Jsc']/jsc_theory) * 100
 
    return {
        'recombination': {'type': recomb_type, 'n': n},
        'series_resistance': {'value': Rs, 'quality': rs_quality},
        'shunt_resistance': {'value': Rsh, 'quality': rsh_quality},
        'losses': {'Voc_loss': voc_loss, 'Jsc_loss': jsc_loss}
    }

class QualityFlags:
    GOOD = "GOOD"
    WARNING = "WARNING"
    FAILED = "FAILED"
 
def validate_physics(self, params):
    issues = []
    if params['Jph'] < params['Jsc_ext']:
        issues.append("Jph < Jsc — 物理矛盾")
    if params['n'] < 1.0:
        issues.append("n < 1.0 — 非物理")
    elif params['n'] > 2.5:
        issues.append("n > 2.5 — 极高度非理想")
    if params['Rs'] > 10:
        issues.append("Rs > 10 — 串联问题严重")
    if params['Rsh'] < 100:
        issues.append("Rsh < 100 — 严重分流")
    return len(issues) == 0, issues

# ===== Step 1: 修改参数文件 =====
def set_params(p):
    # L2_parameters.txt — 4 个层参数
    l2 = os.path.join(WORK_DIR, 'L2_parameters.txt')
    with open(l2, 'r+', encoding='utf-8') as f:
        c = f.read()
        # 正则替换：匹配 "param = 数字" 并替换为新区
        c = re.sub(r'mu_n\s*=\s*[\d.Ee+-]+',      f'mu_n = {p["mu_n"]:.6e}', c)
        c = re.sub(r'mu_p\s*=\s*[\d.Ee+-]+',      f'mu_p = {p["mu_p"]:.6e}', c)
        c = re.sub(r'N_t_bulk\s*=\s*[\d.Ee+-]+',  f'N_t_bulk = {p["Nt_bulk"]:.6e}', c)
        c = re.sub(r'N_t_int\s*=\s*[\d.Ee+-]+',   f'N_t_int = {p["Nt_int"]:.6e}', c)
        f.seek(0); f.write(c); f.truncate()
 
    # simulation_setup.txt — R_shunt
    setup = os.path.join(WORK_DIR, 'simulation_setup.txt')
    with open(setup, 'r+', encoding='utf-8') as f:
        c = f.read()
        c = re.sub(r'R_shunt\s*=\s*-?[\d.Ee+-]+', f'R_shunt = {p["R_shunt"]:.6e}', c)
        f.seek(0); f.write(c); f.truncate()
 
# ===== Step 2: 准备实验 J-V 数据 =====
def prepare_exp_data():
    df = pd.read_csv(os.path.join(KEY_CURVES_DIR, f'{STAGE}_curve.csv'))
    lines = [" Vext Jext"]  # SIMsalabim 要求的文件头
    for _, row in df.iterrows():
        # 电流取负号：实验数据 Jsc<0 (第四象限)，SIMsalabim 用 Jsc>0 (第一象限)
        lines.append(f" {row['Voltage(V)']:12.6E} {-row['Current(mA/cm2)']:12.6E}")
    with open(os.path.join(WORK_DIR, 'JV_Exp.dat'), 'w') as f:
        f.write('\n'.join(lines))
 
# ===== Step 3: 调用 simss.exe =====
def run_simss():
    # 清除上次输出文件，避免误读
    for f in ['JV.dat', 'scPars.dat', 'log.txt', 'Var.dat']:
        fp = os.path.join(WORK_DIR, f)
        if os.path.exists(fp): os.remove(fp)
 
    # 调用 SIMsalabim，90 秒超时
    r = subprocess.run(
        [SIMSS_EXE, 'simulation_setup.txt'],
        cwd=WORK_DIR,
        capture_output=True, text=True,
        timeout=90, encoding='utf-8', errors='replace'
    )
 
    # 检查是否崩溃
    if r.returncode != 0:
        for line in (r.stdout + r.stderr).split('\n'):
            if 'Cannot' in line:
                return None  # SIMsalabim 无法求解此参数组合
 
    # 从 stdout 解析 fitError
    for line in r.stdout.split('\n'):
        if 'fitError:' in line:
            return float(line.split(':')[1].strip().split()[0])
    return None

mu_n_vals = np.logspace(np.log10(1e-6), np.log10(5e-3), 6)
mu_p_vals = np.logspace(np.log10(1e-7), np.log10(5e-3), 6)
Nt_vals   = np.logspace(np.log10(1e18),  np.log10(1e24), 5)
Nt_int_vals = np.logspace(np.log10(1e10), np.log10(1e16), 5)
Rsh_vals  = np.logspace(np.log10(1e-4),  np.log10(1e4),  5)
 
best = {'fitError': float('inf')}
for mu_n in mu_n_vals:
    for mu_p in mu_p_vals:
        for nt in Nt_vals:
            for n_int in Nt_int_vals:
                for rsh in Rsh_vals:
                    set_params({...}); e = run_simss()
                    if e and e < best['fitError']:
                        best = {...}; best['fitError'] = e

def coord_descent_clamped(start, max_cyc=5, dr=1.0, n=15):
    """带物理约束的坐标下降"""
    best = clamp_params(dict(start))      # Step 0: 物理约束
    set_params(best)
    best['fitError'] = run_simss()
 
    for cyc in range(max_cyc):
        improved = False
        for pn in PARAM_NAMES:             # 逐参数优化
            log_val = np.log10(max(best[pn], 1e-30))
 
            # 搜索范围 = 当前值 ±dr 对数空间，并 clamp 到物理边界
            lo_log = max(np.log10(BOUNDS[pn][0]), log_val - dr)
            hi_log = min(np.log10(BOUNDS[pn][1]), log_val + dr)
            if hi_log - lo_log < 0.01: continue  # 搜索空间太窄，跳过
 
            vals = np.linspace(lo_log, hi_log, n)  # 对数均匀取 n 点
            for v in 10**vals:
                trial = dict(best); trial[pn] = v
                trial = clamp_params(trial)         # 每次都要约束！
                set_params(trial)
                e = run_simss()
                if e is not None and e < best['fitError']:
                    best[pn] = v
                    best['fitError'] = e
                    improved = True
 
        if not improved:      # 本轮无任何参数改善 → 已收敛
            break
    return best

# ===== 物理硬约束边界 =====
BOUNDS = {
    'mu_n':    (1e-7, 5e-3),      # max 50 cm²/V·s
    'mu_p':    (1e-7, 5e-3),
    'Nt_bulk': (1e18, 1e24),      # 接近 Nc/2 上限
    'Nt_int':  (1e10, 1e16),
    'R_shunt': (1e-4, 1e4),
}
 
# ===== 迁移率比值约束 =====
RATIO_MIN = 0.001   # μp ≥ μn/1000
RATIO_MAX = 2.0     # μp ≤ 2×μn
 
def clamp_params(p):
    # 1. 边界约束
    for k in PARAM_NAMES:
        lo, hi = BOUNDS[k]
        p[k] = max(lo, min(hi, p[k]))
 
    # 2. 比值约束（防 μp/μn 极端不合理）
    if p['mu_p'] / p['mu_n'] < RATIO_MIN:
        p['mu_p'] = p['mu_n'] * RATIO_MIN
    if p['mu_p'] / p['mu_n'] > RATIO_MAX:
        p['mu_p'] = p['mu_n'] * RATIO_MAX
 
    return p

def random_start():
    """在对数空间中均匀随机采样"""
    p = {}
    for name in PARAM_NAMES:
        lo, hi = BOUNDS[name]
        log_val = np.random.uniform(np.log10(lo), np.log10(hi))
        p[name] = 10**log_val
    return clamp_params(p)
 
# 两阶段优化
def run_optimization(name, start):
    # Phase 1: ±1 decade, 5 cycles — 粗搜索
    p1 = coord_descent_clamped(start, max_cyc=5, dr=1.0, n=15)
    # Phase 2: ±2 decades, 3 cycles — 扩大搜索
    p2 = coord_descent_clamped(p1, max_cyc=3, dr=2.0, n=15)
    return p2 if p2['fitError'] < p1['fitError'] else p1

# 并行运行 5 个阶段
python exp5_constrained.py initial &
python exp5_constrained.py early   &
python exp5_constrained.py mid     &
python exp5_constrained.py late    &
python exp5_constrained.py final   &

WORK_DIR = f"G:\\OpenClaw-Workspace\\simsalabim_parallel\\{STAGE}_exp5"
 
def setup_workspace():
    shutil.rmtree(WORK_DIR)                           # 清除旧副本
    shutil.copytree(ORIGINAL_DIR, WORK_DIR)            # 创建独立副本
    shutil.copytree(data_src, WORK_DIR + "\\Data")     # 复制材料数据

STAGES = ['initial', 'early', 'mid', 'late', 'final']
 
all_results = {}
for stage in STAGES:
    # 1. 加载实验 J-V（来自 key_curves 目录的 CSV）
    df = pd.read_csv(f'{KEY_CURVES_DIR}/{stage}_curve.csv')
 
    # 2. 写入 SIMsalabim 输入文件（格式：Vext Jext）
    lines = [" Vext Jext"]
    for _, row in df.iterrows():
        lines.append(f" {row['Voltage(V)']:12.6E} {-row['Current(mA/cm2)']:12.6E}")
    with open(f'{WORK_DIR}/JV_Exp.dat', 'w') as f:
        f.write('\n'.join(lines))
 
    # 3. 运行坐标下降（exp5 最优方法）
    result = coord_descent_clamped(start_params[stage], max_cyc=3, dr=1.0, n=15)
 
    # 4. 保存各阶段结果
    all_results[stage] = result
    with open(f'results/{stage}_results.json', 'w') as f:
        json.dump(result, f)

# 从两个阶段的结果计算参数变化
def compute_degradation_metrics(initial, final):
    metrics = {
        # 各参数相对变化
        'mu_n_degradation': (initial['mu_n'] - final['mu_n']) / initial['mu_n'],
        'Nt_bulk_increase': final['Nt_bulk'] / initial['Nt_bulk'],
        'Nt_int_increase': final['Nt_int'] / initial['Nt_int'],
        'PCE_loss': (initial['PCE'] - final['PCE']) / initial['PCE'],
        # fitError 对比
        'fitError_improvement': (exp1_fitError - optimized_fitError),
    }
    return metrics

# 读取两个数据集的结果
ds1 = load_results('aging_analysis_r1-r4/results/')
ds2 = load_results('aging_analysis_v2/results/')
 
# 关键对比
comparisons = {
    'PCE_loss':        {'ds1': 0.881, 'ds2': 0.332},
    'mu_n_degradation': {'ds1': 0.96,  'ds2': 0.57},
    'Nt_int_increase':  {'ds1': 50,    'ds2': 2.2},
    'n_final':          {'ds1': 2.5,   'ds2': 3.0},
}