#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Hall効果移動度解析・フィッティングツール: tkminfit callback(q) 版

主な変更点:
- scipy.optimize.minimize の直呼びをやめ、tkminfit.minimize_lsq を使用。
- callback は tkminfit の新仕様 callback(q) に対応。
- callback 内の iteration 管理、unpack(q)、objective(q)、リアルタイム描画、GIF用frame保存はアプリ側で実装。
- optid=1 は optid_lin の値に関係なくすべて q に含める。
- 元コード同様、y=log10(mu) とし、log10移動度残差を最小化。
"""

from __future__ import annotations

import argparse
import json
import os
import sys
import traceback
from pathlib import Path
from typing import Dict, Mapping, Sequence, Tuple

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter

# 必要ならローカル開発環境の tklsq パスを追加
_tkprog_path = os.environ.get("tkprog_X_path")
if _tkprog_path:
    sys.path.append(os.path.join(_tkprog_path, "regression"))

import tklsq.tkdataio as tkdataio
import tklsq.tkparamio as tkparamio
import tklsq.tklsq as tklsq_core
import tklsq.tknlsq as tknlsq
import tklsq.tkfitdiag as tkfitdiag
import tklsq.tkminfit as tkminfit


K_B = 8.617333262e-5  # eV/K
PARAM_NAMES = ["aop", "a1", "a2", "a3", "VB"]

DEFAULT_PARAMS = [
    {"varname": "aop", "optid": 1, "optid_lin": 0, "p0": 1e-3, "pscale": "log", "dp": 1e-6, "pmin": 1e-300, "pmax": "", "kpenalty": 1e12},
    {"varname": "a1",  "optid": 1, "optid_lin": 0, "p0": 1e-7, "pscale": "log", "dp": 1e-10, "pmin": 1e-300, "pmax": "", "kpenalty": 1e12},
    {"varname": "a2",  "optid": 1, "optid_lin": 0, "p0": 1e-4, "pscale": "log", "dp": 1e-7, "pmin": 1e-300, "pmax": "", "kpenalty": 1e12},
    {"varname": "a3",  "optid": 1, "optid_lin": 0, "p0": 1e1,  "pscale": "log", "dp": 1e-2, "pmin": 1e-300, "pmax": "", "kpenalty": 1e12},
    {"varname": "VB",  "optid": 1, "optid_lin": 0, "p0": 0.0,  "pscale": "linear", "dp": 1e-4, "pmin": "", "pmax": "", "kpenalty": 1e12},
]


def stem_from_infile(infile: str | Path) -> str:
    return str(Path(infile).with_suffix(""))


def auto_paramfile(infile: str | Path) -> str:
    return stem_from_infile(infile) + "_params.csv"


def auto_outfile(infile: str | Path) -> str:
    return stem_from_infile(infile)


def save_log(path: str | Path) -> Path:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    print(f"[save] writing {path}")
    return path


def save_done(path: str | Path) -> None:
    print(f"[save] done {Path(path)}")


def read_xy_hall(args) -> Tuple[np.ndarray, np.ndarray, str, str]:
    x, y, xlabel, ylabel = tkdataio.read_xy(
        args.infile,
        x=args.temp_col,
        y=args.mu_col,
        sheet_name=args.sheet_name,
        xmin=args.Tfitmin_read,
        xmax=args.Tfitmax_read,
        dropna=True,
    )
    mask = np.isfinite(x) & np.isfinite(y) & (x > 0) & (y > 0)
    if not np.all(mask):
        print(f"[read] dropped {np.size(mask) - int(np.sum(mask))} invalid rows")
    return x[mask], y[mask], xlabel, ylabel


def get_inv_mu_components(T: np.ndarray, p: Mapping[str, float], Eop: float):
    T = np.asarray(T, dtype=float)
    aop = float(p["aop"])
    a1 = float(p["a1"])
    a2 = float(p["a2"])
    a3 = float(p["a3"])
    VB = float(p["VB"])

    f_op = 1.0 / (np.exp(np.clip(Eop / (K_B * T), -700.0, 700.0)) - 1.0)
    f_ac = T ** 1.5
    f_ni = np.ones_like(T)
    f_ii = T ** -1.5

    components = {
        "Optical Phonon": aop * f_op,
        "Acoustic Phonon": a1 * f_ac,
        "Neutral Impurity": a2 * f_ni,
        "Ionized Impurity": a3 * f_ii,
    }

    inv_mu_bulk = np.maximum(sum(components.values()), 1e-300)
    expo = np.clip(VB / (K_B * T), -700.0, 700.0)
    inv_mu_total = inv_mu_bulk * np.exp(expo)
    components["Grain Boundary"] = inv_mu_total - inv_mu_bulk

    return components, np.maximum(inv_mu_total, 1e-300)


def mobility_model(T: np.ndarray, p: Mapping[str, float]) -> np.ndarray:
    _, inv_total = get_inv_mu_components(T, p, float(p.get("Eop", 0.045)))
    return 1.0 / np.maximum(inv_total, 1e-300)


def model(x: np.ndarray, p: Mapping[str, float]) -> np.ndarray:
    """fit用モデル。元コードに合わせて log10(mu_model) を返す。"""
    mu = mobility_model(x, p)
    return np.log10(np.maximum(mu, 1e-300))


def lfit_design_matrix(T: np.ndarray, Eop: float) -> np.ndarray:
    T = np.asarray(T, dtype=float)
    f_op = 1.0 / (np.exp(np.clip(Eop / (K_B * T), -700.0, 700.0)) - 1.0)
    f_ac = T ** 1.5
    f_ni = np.ones_like(T)
    f_ii = T ** -1.5
    return np.column_stack([f_op, f_ac, f_ni, f_ii])


def fit_mask(T: np.ndarray, tmin: float, tmax: float) -> np.ndarray:
    return (T >= tmin) & (T <= tmax)


def values_from_params(params) -> Dict[str, float]:
    return {name: float(row["p0"]) for name, row in params.items()}


def optid1_names(params) -> list[str]:
    """mode=fit用。optid=1 は optid_lin の値に関係なくすべて q に含める。"""
    return [name for name, row in params.items() if int(row.get("optid", 0)) == 1]


def pscale_of(params, name: str) -> str:
    return str(params[name].get("pscale", "linear")).strip().lower() or "linear"

def positive_floor_from_param(params, name: str, default: float = 1e-300) -> float:
    """
    pscale=log 用の最小正値を返す。

    優先順位:
        1. pmin (>0)
        2. default
    """
    row = params[name]

    try:
        pmin = float(row.get("pmin", np.nan))
        if np.isfinite(pmin) and pmin > 0:
            return float(pmin)
    except Exception:
        pass

    return float(default)
    
def pack(values: Mapping[str, float], names: Sequence[str], params) -> np.ndarray:
    q = []

    for name in names:
        v = float(values[name])

        if pscale_of(params, name) == "log":

            floor = positive_floor_from_param(params, name)

            if not np.isfinite(v) or v <= 0:
                print(
                    f"[warn] {name}={v:g} is invalid for log scale; "
                    f"clipped to {floor:g}"
                )
                v = floor

            v = max(v, floor)

            q.append(np.log10(v))

        else:
            q.append(v)

    return np.asarray(q, dtype=float)


def unpack_factory(base_values: Mapping[str, float], names: Sequence[str], params, eop: float):
    def unpack(q: Sequence[float]) -> Dict[str, float]:
        p = {k: float(v) for k, v in base_values.items()}
        for name, qv in zip(names, np.asarray(q, dtype=float).reshape(-1)):
            if pscale_of(params, name) == "log":
                p[name] = float(10.0 ** qv)
            else:
                p[name] = float(qv)
        p["Eop"] = float(eop)
        return p

    return unpack


def apply_cli_defaults(defaults, args):
    rows = [dict(r) for r in defaults]
    fixed = [s.strip() for s in str(args.fix).split(",") if s.strip()]
    for row in rows:
        name = row["varname"]
        if name in fixed:
            row["optid"] = 0
        cli_attr = f"p0_{name}"
        if hasattr(args, cli_attr) and getattr(args, cli_attr) is not None:
            row["p0"] = getattr(args, cli_attr)
    return rows


def read_params(args):
    defaults = apply_cli_defaults(DEFAULT_PARAMS, args)
    path = args.paramfile or auto_paramfile(args.infile)
    created = not Path(path).exists()
    if created:
        print(f"[save] writing {path}")
    params = tkparamio.read_param_csv(path, defaults=defaults, create_if_missing=True)
    if created:
        print(f"[save] done {path}")
    return params, path


def prediction_band_delta(
    x: np.ndarray,
    p_fit: Mapping[str, float],
    free_names: Sequence[str],
    cov,
    params,
    nsigma: float,
    rel_step: float,
    abs_step: float,
):
    y0 = model(x, p_fit)
    if cov is None or len(free_names) == 0:
        return y0, None, None, None

    q0 = pack(p_fit, free_names, params)
    unpack = unpack_factory(p_fit, free_names, params, p_fit.get("Eop", 0.045))

    def y_of_q(q):
        return model(x, unpack(q))

    G = np.zeros((len(x), len(free_names)), dtype=float)
    for j in range(len(free_names)):
        dq = rel_step * (abs(q0[j]) + 1.0) + abs_step
        qp = q0.copy()
        qm = q0.copy()
        qp[j] += dq
        qm[j] -= dq
        G[:, j] = (y_of_q(qp) - y_of_q(qm)) / (2.0 * dq)

    var_y = np.einsum("ni,ij,nj->n", G, cov, G)
    sigma_y = np.sqrt(np.maximum(var_y, 0.0))
    return y0, y0 - nsigma * sigma_y, y0 + nsigma * sigma_y, sigma_y


def setup_fit_axes(title: str):
    fig, (ax1, ax2) = plt.subplots(
        2,
        1,
        figsize=(8, 7),
        sharex=True,
        gridspec_kw={"height_ratios": [3, 1]},
    )
    ax1.set_title(title)
    ax1.set_ylabel("log10 Mobility")
    ax1.grid(True, alpha=0.3)
    ax2.set_xlabel("Temperature (K)")
    ax2.set_ylabel("residual")
    ax2.grid(True, alpha=0.3)
    return fig, ax1, ax2


def draw_fit_residual(fig, ax1, ax2, x, y, y_fit, residual, title=None, y_lower=None, y_upper=None):
    ax1.clear()
    ax2.clear()
    idx = np.argsort(x)
    ax1.scatter(x, y, s=24, color="tab:red", alpha=0.7, label="data")
    ax1.plot(x[idx], y_fit[idx], color="tab:blue", lw=2, label="fit")
    if y_lower is not None and y_upper is not None:
        ax1.fill_between(x[idx], y_lower[idx], y_upper[idx], color="tab:blue", alpha=0.18, label="±1σ")
    ax1.set_ylabel("log10 Mobility")
    if title:
        ax1.set_title(title)
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    ax2.axhline(0.0, color="black", lw=1)
    ax2.scatter(x, residual, s=22, color="tab:purple", alpha=0.75)
    ax2.set_xlabel("Temperature (K)")
    ax2.set_ylabel("residual")
    ax2.grid(True, alpha=0.3)
    fig.tight_layout()


def save_fit_png(path, x, y, y_fit, residual, title, show, y_lower=None, y_upper=None):
    fig, ax1, ax2 = setup_fit_axes(title)
    draw_fit_residual(fig, ax1, ax2, x, y, y_fit, residual, title, y_lower, y_upper)
    p = save_log(path)
    fig.savefig(p, dpi=160)
    save_done(p)
    if int(show) == 1:
        plt.show()
    plt.close(fig)


def save_animation(path, x, y, frames):
    if not frames:
        print("[animation] no frames; skipped")
        return

    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    print(f"[animation] saving animation to {path}")
    print(f"[save] writing {path}")

    fig, ax1, ax2 = setup_fit_axes("fit progress")

    def update(i):
        fr = frames[i]
        draw_fit_residual(
            fig,
            ax1,
            ax2,
            x,
            y,
            fr["y_fit"],
            fr["residual"],
            title=f"iteration={fr['iteration']} objective={fr['objective']:.6g}",
        )
        return []

    ani = FuncAnimation(fig, update, frames=len(frames), interval=300, blit=False)
    ani.save(path, writer=PillowWriter(fps=3))
    plt.close(fig)

    print(f"[save] done {path}")
    print(f"[animation] saved animation to {path}")


def write_excel(path, tables):
    p = save_log(path)
    tkdataio.write_excel_tables(p, tables, index=False)
    save_done(p)


def write_params(path, params, values=None, stderr=None):
    p = save_log(path)
    tkparamio.write_param_csv(p, params, values=values, stderr=stderr)
    save_done(p)


def mode_read(args, x, mu):
    y = np.log10(np.maximum(mu, 1e-300))
    save_fit_png(
        f"{args.outfile}_read.png",
        x,
        y,
        y,
        np.zeros_like(y),
        "Experimental Data",
        args.show,
    )


def mode_lfit(args, x, mu, params, paramfile):
    Eop = float(args.eop)
    X = lfit_design_matrix(x, Eop)
    y_inv = 1.0 / np.maximum(mu, 1e-300)
    res = tklsq_core.linear_lsq(X, y_inv)
    coeffs = np.asarray(res.beta, dtype=float)

    for i, name in enumerate(["aop", "a1", "a2", "a3"]):
        if pscale_of(params, name) == "log":
            floor = positive_floor_from_param(params, name)
            if not np.isfinite(coeffs[i]) or coeffs[i] <= 0:
                print(
                    f"[lfit] {name}={coeffs[i]:g} "
                    f"clipped to {floor:g}"
                )
                coeffs[i] = floor

            coeffs[i] = max(coeffs[i], floor)
        else:
            if not np.isfinite(coeffs[i]):
                coeffs[i] = 0.0

    values = values_from_params(params)
    for name, v in zip(["aop", "a1", "a2", "a3"], coeffs):
        values[name] = float(v)
    values["VB"] = float(values.get("VB", 0.0))
    values["Eop"] = Eop

    stderr = {name: None for name in params.keys()}
    if res.beta_std is not None:
        for name, se in zip(["aop", "a1", "a2", "a3"], res.beta_std):
            stderr[name] = float(se)

    print("[lfit] beta=", coeffs)
    if res.warning:
        print(res.warning)

    if int(args.save) == 1:
        write_params(paramfile, params, values=values, stderr=stderr)

    y = np.log10(np.maximum(mu, 1e-300))
    y_fit = model(x, values)
    if int(args.save) == 1:
        save_fit_png(
            f"{args.outfile}_lfit.png",
            x,
            y,
            y_fit,
            y - y_fit,
            "LLSQ Initial Fit (VB fixed)",
            args.show,
        )


def make_stderr_value_space(stderr_q, p_fit, free_names, params):
    """tkminfitのstderrはq空間なので、log scaleは値空間の近似stderrへ変換する。"""
    stderr = {name: None for name in params.keys()}
    if stderr_q is None:
        return stderr

    for name, se in zip(free_names, stderr_q):
        if se is None:
            stderr[name] = None
        elif pscale_of(params, name) == "log":
            stderr[name] = float(abs(p_fit[name]) * np.log(10.0) * se)
        else:
            stderr[name] = float(se)
    return stderr


def mode_fit(args, x_all, mu_all, params, paramfile):
    mfit = fit_mask(x_all, args.Tfitmin, args.Tfitmax)
    x = x_all[mfit]
    mu = mu_all[mfit]
    if x.size < 3:
        raise ValueError("フィットに使える点が少なすぎます。")

    y = np.log10(np.maximum(mu, 1e-300))
    base_values = values_from_params(params)
    base_values["Eop"] = float(args.eop)

    # optid=1 は optid_lin に関係なくすべて q に含める。
    free_names = optid1_names(params)
    unpack = unpack_factory(base_values, free_names, params, float(args.eop))
    q0 = pack(base_values, free_names, params)

    options = {"maxiter": int(args.maxiter)}

    fig, ax1, ax2 = setup_fit_axes("real-time fit")

    def objective(q):
        p = unpack(q)
        r = y - model(x, p)
        penalty = tkparamio.bounds_penalty(params, p)
        return float(r @ r + penalty)

    def residual_func_from_q_space(p_dummy):
        """
        tkminfit.minimize_lsq は通常 p dict を受けるが、ここでは
        free_names と初期値を q 空間に差し替えるため、p_dummy の値を
        q として解釈する。

        callback(q) との一貫性を保つため、tkminfit 内部の q は
        すでに pack 後の内部最適化ベクトルにする。
        """
        q = np.array([float(p_dummy[n]) for n in free_names], dtype=float)
        p = unpack(q)
        return y - model(x, p)

    # tkminfit に渡す params は q 空間用に変換する。
    # 注意:
    #   tkminfit の callback(q) に渡る q は、この qparams の p0 から始まる
    #   内部最適化ベクトルであり、物理パラメータ dict ではない。
    qparams = {}
    for name, qv in zip(free_names, q0):
        row = dict(params[name])
        row["p0"] = float(qv)
        row["pscale"] = "linear"
        row["pmin"] = ""
        row["pmax"] = ""
        qparams[name] = row

    iteration_state = {"i": 0, "frames": []}

    def callback(q):
        iteration_state["i"] += 1
        p = unpack(q)
        obj = objective(q)

        print(f"[callback] iteration={iteration_state['i']} objective={obj:.10g}")
        print("  " + ", ".join(f"{k}={v:.8g}" for k, v in p.items() if k != "Eop"))

        y_fit = model(x, p)
        residual = y - y_fit
        draw_fit_residual(
            fig,
            ax1,
            ax2,
            x,
            y,
            y_fit,
            residual,
            title=f"iteration={iteration_state['i']} objective={obj:.6g}",
        )
        plt.pause(0.01)

        if iteration_state["i"] % args.nplot_interval == 0:
            iteration_state["frames"].append({
                "iteration": iteration_state["i"],
                "params": dict(p),
                "objective": obj,
                "y_fit": y_fit.copy(),
                "residual": residual.copy(),
            })

    res = tkminfit.minimize_lsq(
        residual_func_from_q_space,
        qparams,
        free_names=free_names,
        method=args.method,
        use_penalty=False,
        callback=callback,
        options=options,
        rel_step=args.jac_relstep,
        abs_step=args.jac_absstep,
    )
    plt.close(fig)

    q_fit = np.asarray(res.params_free, dtype=float).reshape(-1)
    p_fit = unpack(q_fit)
    y_fit = model(x, p_fit)
    residual = y - y_fit
    rss = float(residual @ residual)

    print(f"[fit] success={res.success} message={res.message}")
    print(f"[fit] RSS={rss:.10g}")

    # tkminfit が q 空間で共分散を推定する。
    J = res.jacobian
    cov = res.cov_free
    sigma2 = res.sigma2_resid
    dof = res.dof
    stderr = make_stderr_value_space(res.stderr_free, p_fit, free_names, params)

    if res.warning:
        print(res.warning)

    if cov is not None:
        diag = tkfitdiag.diagnose_covariance(
            free_names,
            [q_fit[i] for i in range(len(free_names))],
            cov,
            jacobian=J,
        )
        print("[diagnostics] cond(JTJ)=", diag.cond_jtj)
        print("[diagnostics] eig(JTJ)=", diag.eig_jtj_values_asc)
        print("[diagnostics] corr=")
        print(diag.corr)
        suggestions = tkfitdiag.propose_fix_candidates_from_diagnostics(diag)
    else:
        diag = None
        suggestions = []

    y0, ylo, yhi, sigy = prediction_band_delta(
        x,
        p_fit,
        free_names,
        cov,
        params,
        args.band_sigma,
        args.jac_relstep,
        args.jac_absstep,
    )

    mu_fit = 10.0 ** y0
    mu_lower = 10.0 ** ylo if ylo is not None else np.full_like(mu_fit, np.nan)
    mu_upper = 10.0 ** yhi if yhi is not None else np.full_like(mu_fit, np.nan)

    if int(args.save) == 1:
        write_params(paramfile, params, values=p_fit, stderr=stderr)
        save_fit_png(
            f"{args.outfile}_fit.png",
            x,
            y,
            y0,
            y - y0,
            f"Final Fit (±{args.band_sigma}σ)",
            args.show,
            ylo,
            yhi,
        )
        save_animation(f"{args.outfile}_fit_animation.gif", x, y, iteration_state["frames"])

        result_df = pd.DataFrame({
            "T": x,
            "mu_exp": mu,
            "log10_mu_exp": y,
            "y_fit": mu_fit,
            "y_lower": mu_lower,
            "y_upper": mu_upper,
            "residual": y - y0,
            "log10_y_fit": y0,
            "log10_y_lower": ylo if ylo is not None else np.nan,
            "log10_y_upper": yhi if yhi is not None else np.nan,
        })

        params_df = pd.DataFrame({
            "varname": list(params.keys()),
            "value": [p_fit.get(k, np.nan) for k in params.keys()],
            "stderr": [stderr.get(k, np.nan) if stderr is not None else np.nan for k in params.keys()],
        })

        diag_tables = {"fit": result_df, "params": params_df}
        if diag is not None:
            diag_tables["corr"] = pd.DataFrame(diag.corr, index=free_names, columns=free_names).reset_index()
            diag_tables["eig_jtj"] = pd.DataFrame({"eig_jtj": diag.eig_jtj_values_asc})
            diag_tables["suggestions"] = pd.DataFrame([
                s.__dict__ if hasattr(s, "__dict__") else s
                for s in suggestions
            ])
        write_excel(f"{args.outfile}_fit.xlsx", diag_tables)

        json_path = save_log(f"{args.outfile}_fit_diagnostics.json")
        payload = {
            "success": bool(res.success),
            "message": str(res.message),
            "RSS": rss,
            "sigma2": sigma2,
            "dof": dof,
            "free_names": free_names,
            "params": {k: float(v) for k, v in p_fit.items() if k != "Eop"},
            "q_fit": q_fit.tolist(),
            "cond_JTJ": None if diag is None else diag.cond_jtj,
            "eig_JTJ": None if diag is None or diag.eig_jtj_values_asc is None else diag.eig_jtj_values_asc.tolist(),
        }
        with open(json_path, "w", encoding="utf-8") as f:
            json.dump(payload, f, indent=2, ensure_ascii=False)
        save_done(json_path)


def mode_weight(args, x, mu, params):
    p = values_from_params(params)
    p["Eop"] = float(args.eop)
    components, total = get_inv_mu_components(x, p, float(args.eop))

    fig, ax = plt.subplots(figsize=(9, 6))
    idx = np.argsort(x)
    for name, inv_mu in components.items():
        weight = 100.0 * inv_mu / np.maximum(total, 1e-300)
        ax.plot(x[idx], weight[idx], marker="o", ms=4, label=name)

    ax.set_xlabel("Temperature (K)")
    ax.set_ylabel("Contribution to Scattering (%)")
    ax.set_ylim(-5, 105)
    ax.grid(True, alpha=0.3)
    ax.legend()
    fig.tight_layout()

    if int(args.save) == 1:
        p_png = save_log(f"{args.outfile}_weight.png")
        fig.savefig(p_png, dpi=160)
        save_done(p_png)
    if int(args.show) == 1:
        plt.show()
    plt.close(fig)


def mode_sim(args, x, mu, params):
    p = values_from_params(params)
    p["Eop"] = float(args.eop)
    y_data = np.log10(np.maximum(mu, 1e-300))
    y_fit = model(x, p)
    if int(args.save) == 1:
        save_fit_png(
            f"{args.outfile}_sim.png",
            x,
            y_data,
            y_fit,
            y_data - y_fit,
            "Simulation from params",
            args.show,
        )


def mode_model_select(args, x, mu, params):
    print("[model_select] このHall移動度モデルでは、現在は全パラメータモデルのみを評価します。")
    mode_fit(args, x, mu, params, args.paramfile or auto_paramfile(args.infile))


def build_parser():
    parser = argparse.ArgumentParser(description="Hall効果移動度フィッティング tkminfit callback(q) 版")
    parser.add_argument("--mode", required=True, choices=["read", "sim", "lfit", "llsq", "fit", "model_select", "weight"])
    parser.add_argument("--infile", required=True, help="入力ファイル CSV/TSV/TXT/XLSX")
    parser.add_argument("--save", type=int, choices=[0, 1], default=1)
    parser.add_argument("--show", type=int, choices=[0, 1], default=1)
    parser.add_argument("--paramfile", default=None)
    parser.add_argument("--outfile", default=None)
    parser.add_argument("--nplot_interval", type=int, default=10)
    parser.add_argument("--method", default="Nelder-Mead")
    parser.add_argument("--maxiter", type=int, default=5000)
    parser.add_argument("--temp_col", default=0, help="温度列 index または列名")
    parser.add_argument("--mu_col", default=2, help="移動度列 index または列名")
    parser.add_argument("--sheet_name", default=0)
    parser.add_argument("--eop", type=float, default=0.045, help="光学フォノンエネルギー eV")
    parser.add_argument("--fix", default="", help="CSV未作成時のみ有効。固定するパラメータ名をカンマ区切り指定")
    parser.add_argument("--Tfitmin", type=float, default=-1e100)
    parser.add_argument("--Tfitmax", type=float, default=+1e100)
    parser.add_argument("--Tfitmin_read", type=float, default=-1e100)
    parser.add_argument("--Tfitmax_read", type=float, default=+1e100)
    parser.add_argument("--band_sigma", type=float, default=1.0)
    parser.add_argument("--jac_relstep", type=float, default=1e-6)
    parser.add_argument("--jac_absstep", type=float, default=1e-12)
    parser.add_argument("--p0_aop", type=float, default=None)
    parser.add_argument("--p0_a1", type=float, default=None)
    parser.add_argument("--p0_a2", type=float, default=None)
    parser.add_argument("--p0_a3", type=float, default=None)
    parser.add_argument("--p0_VB", type=float, default=None)
    return parser


def main():
    try:
        parser = build_parser()
        args = parser.parse_args()

        if args.outfile is None:
            args.outfile = auto_outfile(args.infile)
        if args.paramfile is None:
            args.paramfile = auto_paramfile(args.infile)

        print(f"infile={args.infile}")
        print(f"outfile={args.outfile}")
        print(f"paramfile={args.paramfile}")
        print(f"mode={args.mode}")
        print(f"method={args.method}")
        print(f"eop={args.eop} eV")

        x, mu, xlabel, ylabel = read_xy_hall(args)
        params, paramfile = read_params(args)

        mode = "lfit" if args.mode == "llsq" else args.mode
        if mode == "read":
            mode_read(args, x, mu)
        elif mode == "lfit":
            mode_lfit(args, x, mu, params, paramfile)
        elif mode == "fit":
            mode_fit(args, x, mu, params, paramfile)
        elif mode == "sim":
            mode_sim(args, x, mu, params)
        elif mode == "weight":
            mode_weight(args, x, mu, params)
        elif mode == "model_select":
            mode_model_select(args, x, mu, params)
        else:
            raise ValueError(f"unknown mode: {args.mode}")
    except Exception:
        traceback.print_exc()
        raise


if __name__ == "__main__":
    main()