ベイズ最適化による材料探索｜Pythonで次の実験条件を選ぶ方法

"""Educational, finite-candidate, single-objective Bayesian optimization.
Python 3.12; numpy, pandas, scipy, scikit-learn, matplotlib.
All demo data are synthetic, not recommended ceramic processing conditions.
"""
import argparse
import json
import warnings
from pathlib import Path
import numpy as np
import pandas as pd
import scipy
from scipy.stats import norm
import sklearn
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, Matern
from sklearn.exceptions import ConvergenceWarning
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

FEATURES = ['temperature_c', 'hold_time_min']
TARGET = 'strength_mpa'
LOWER = np.array([1000., 30.])
UPPER = np.array([1400., 150.])
SEED = 42

def acquisition(mean, std, best, xi=0.0, kappa=2.0):
    mean, std = np.asarray(mean), np.asarray(std)
    improvement = mean - best - xi
    positive = std > 1e-12
    z = np.divide(improvement, std, out=np.zeros_like(mean), where=positive)
    ei = np.where(positive, improvement * norm.cdf(z) + std * norm.pdf(z),
                  np.maximum(improvement, 0.0))
    pi = np.where(positive, norm.cdf(z), (improvement > 0).astype(float))
    return np.maximum(ei, 0.0), pi, mean + kappa * std

def validate(observations, candidates):
    if not {'experiment_id', TARGET, *FEATURES}.issubset(observations.columns):
        raise ValueError('Observations need experiment_id, strength_mpa and both feature columns.')
    if not {'candidate_id', *FEATURES}.issubset(candidates.columns):
        raise ValueError('Candidates need candidate_id and both feature columns.')
    observations, candidates = observations.copy(), candidates.copy()
    for frame, id_col in [(observations, 'experiment_id'), (candidates, 'candidate_id')]:
        if frame[id_col].isna().any() or frame[id_col].astype(str).str.strip().eq('').any():
            raise ValueError('Empty IDs are not allowed.')
        if frame[id_col].astype(str).duplicated().any():
            raise ValueError('IDs must be unique within each CSV.')
        for col in FEATURES + ([TARGET] if id_col == 'experiment_id' else []):
            frame[col] = pd.to_numeric(frame[col], errors='raise')
        values = frame[FEATURES].to_numpy(float)
        if not np.isfinite(values).all():
            raise ValueError('Feature values must be finite; no automatic imputation.')
        if ((values < LOWER) | (values > UPPER)).any():
            raise ValueError('A condition is outside the declared search bounds.')
    if not np.isfinite(observations[TARGET]).all():
        raise ValueError('Measured targets must be finite; failed experiments need separate handling.')
    if len(observations) < 4 or len(observations[FEATURES].drop_duplicates()) < 4:
        raise ValueError('This teaching script requires at least four distinct measured conditions.')
    if observations[FEATURES].nunique().min() < 2:
        raise ValueError('Both variables need variation in the initial data.')
    if candidates[FEATURES].duplicated().any():
        raise ValueError('Duplicate candidate conditions; canonicalize values using instrument resolution.')
    # Exact matching requires canonical values/units (e.g., instrument ticks).
    observed_keys = set(map(tuple, observations[FEATURES].to_numpy(float)))
    keep = [tuple(row) not in observed_keys for row in candidates[FEATURES].to_numpy(float)]
    candidates = candidates.loc[keep].reset_index(drop=True)
    if candidates.empty:
        raise ValueError('No unmeasured candidates remain.')
    return observations, candidates

def propose(observations, candidates, *, minimize=False, noise_std=0.5,
            method='ei', xi=0.0, kappa=2.0):
    if not np.isfinite([noise_std, xi, kappa]).all() or noise_std <= 0 or xi < 0 or kappa < 0:
        raise ValueError('noise_std > 0, xi >= 0 and kappa >= 0 are required.')
    observations, candidates = validate(observations, candidates)
    x = (observations[FEATURES].to_numpy(float) - LOWER) / (UPPER - LOWER)
    xc = (candidates[FEATURES].to_numpy(float) - LOWER) / (UPPER - LOWER)
    sign = -1.0 if minimize else 1.0
    utility = sign * observations[TARGET].to_numpy(float)
    center, scale = utility.mean(), utility.std(ddof=0)
    if scale <= 1e-12:
        raise ValueError('Constant targets: review measurements and initial design before optimization.')
    kernel = ConstantKernel(1.0, (0.01, 100.0)) * Matern(
        length_scale=[0.3, 0.3], length_scale_bounds=(0.03, 3.0), nu=2.5)
    model = GaussianProcessRegressor(kernel=kernel, alpha=(noise_std / scale)**2,
        normalize_y=False, n_restarts_optimizer=2, random_state=SEED)
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter('always', ConvergenceWarning)
        model.fit(x, (utility-center)/scale)
    m, s = model.predict(xc, return_std=True)
    m, s = m*scale+center, s*scale
    ei, pi, ucb = acquisition(m, s, utility.max(), xi, kappa)
    scores = {'ei': ei, 'pi': pi, 'ucb': ucb}[method]
    if not np.isfinite(np.column_stack([m, s, ei, pi, ucb])).all():
        raise ValueError('Non-finite predictions/acquisition values; do not select a candidate.')
    if method == 'ei' and scores.max() <= 0:
        raise ValueError('All EI scores are zero: review scaling/model or use a robust LogEI implementation.')
    result = candidates.copy()
    result['predicted_mean'] = sign*m
    result['latent_std'] = s
    result['ei'] = ei
    result['pi'] = pi
    result['ucb_utility'] = ucb
    result['selection_score'] = scores
    result = result.sort_values('selection_score', ascending=False, kind='stable').reset_index(drop=True)
    result.insert(0, 'rank', np.arange(1, len(result)+1))
    info = {'kernel': str(model.kernel_), 'noise_std_target_units': noise_std,
        'direction': 'minimize' if minimize else 'maximize', 'method': method,
        'xi_target_units': xi, 'kappa': kappa, 'seed': SEED,
        'bounds': {'lower': LOWER.tolist(), 'upper': UPPER.tolist()},
        'warnings': [str(w.message) for w in caught],
        'note': 'latent_std excludes independent observation noise; demo does not validate real processes.'}
    return result, info

def synthetic_truth(frame):
    t, h = frame[FEATURES].to_numpy(float).T
    return 70 + 30*np.exp(-((t-1260)/100)**2-((h-85)/35)**2) - 0.00005*(t-1200)**2

def make_demo():
    candidates = pd.DataFrame([(t,h) for t in range(1000,1401,20)
        for h in range(30,151,10)], columns=FEATURES)
    candidates.insert(0, 'candidate_id', [f'C{i:03}' for i in range(len(candidates))])
    initial = np.array([[1020,40],[1380,40],[1020,140],[1380,140],
                        [1200,90],[1100,90],[1300,90],[1200,50]], float)
    observations = pd.DataFrame(initial, columns=FEATURES)
    observations.insert(0, 'experiment_id', [f'E{i:03}' for i in range(len(initial))])
    observations[TARGET] = synthetic_truth(observations) + np.random.default_rng(SEED).normal(0,0.5,len(initial))
    return observations, candidates

def save_plots(observations, ranked, history, destination):
    fig, ax = plt.subplots(figsize=(7,5))
    pts = ax.scatter(ranked[FEATURES[0]], ranked[FEATURES[1]], c=ranked['ei'], cmap='viridis', s=35)
    ax.scatter(observations[FEATURES[0]], observations[FEATURES[1]], c='white', edgecolors='black', s=65, label='Measured')
    top = ranked.iloc[0]
    ax.scatter(top[FEATURES[0]], top[FEATURES[1]], c='red', marker='*', s=200, label='Next candidate')
    ax.set(xlabel='Temperature (synthetic, C)', ylabel='Hold time (synthetic, min)', title='Finite-candidate EI: synthetic example')
    ax.legend(); fig.colorbar(pts, ax=ax, label='Expected improvement (MPa)')
    fig.tight_layout(); fig.savefig(destination/'candidate_map.png', dpi=150); plt.close(fig)
    if history:
        hist = pd.DataFrame(history)
        hist.to_csv(destination/'simulation_history.csv', index=False)
        fig, ax = plt.subplots(figsize=(7,4))
        ax.plot(hist['n_observed'], hist['best_observed'], 'o-')
        ax.set(xlabel='Number of observations', ylabel='Best observed strength (synthetic, MPa)',
               title='One synthetic run; no guaranteed improvement')
        fig.tight_layout(); fig.savefig(destination/'progress.png', dpi=150); plt.close(fig)

def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--observations'); parser.add_argument('--candidates')
    parser.add_argument('--output', default='bo_results')
    parser.add_argument('--method', choices=['ei','pi','ucb'], default='ei')
    parser.add_argument('--minimize', action='store_true')
    parser.add_argument('--noise-std', type=float, default=0.5)
    parser.add_argument('--xi', type=float, default=0.0)
    parser.add_argument('--kappa', type=float, default=2.0)
    parser.add_argument('--simulate', type=int, default=0, help='Demo-only number of additional synthetic experiments')
    args = parser.parse_args()
    if bool(args.observations) != bool(args.candidates):
        parser.error('Supply both --observations and --candidates, or neither for the synthetic demo.')
    if args.simulate < 0 or (args.simulate and args.observations):
        parser.error('--simulate is non-negative and available only for the built-in synthetic demo.')
    observations, candidates = (pd.read_csv(args.observations), pd.read_csv(args.candidates)) if args.observations else make_demo()
    destination = Path(args.output); destination.mkdir(parents=True, exist_ok=True)
    observations.to_csv(destination/'observations_input.csv', index=False)
    candidates.to_csv(destination/'candidates_input.csv', index=False)
    settings = dict(minimize=args.minimize, noise_std=args.noise_std, method=args.method, xi=args.xi, kappa=args.kappa)
    first_ranked, first_info = propose(observations, candidates, **settings)
    first_ranked.to_csv(destination/'ranked_candidates_initial.csv', index=False)
    initial_observations = observations.copy()
    best = lambda f: float(f[TARGET].min() if args.minimize else f[TARGET].max())
    history = [{'n_observed':len(observations), 'best_observed':best(observations)}]
    rng = np.random.default_rng(SEED+1)
    for step in range(args.simulate):
        ranked, info = propose(observations, candidates, **settings)
        row = ranked.iloc[0]
        added = pd.DataFrame([{ 'experiment_id':f'SIM{step:03}', **{c:row[c] for c in FEATURES}}])
        added[TARGET] = synthetic_truth(added) + rng.normal(0,0.5,len(added))
        observations = pd.concat([observations, added], ignore_index=True)
        history.append({'n_observed':len(observations), 'best_observed':best(observations)})
    ranked, info = propose(observations, candidates, **settings)
    ranked.to_csv(destination/'ranked_candidates.csv', index=False)
    ranked.head(1).to_csv(destination/'next_candidate.csv', index=False)
    observations.to_csv(destination/'observations_updated.csv', index=False)
    info.update({'python_note':'Tested environment recorded in article', 'numpy':np.__version__,
        'pandas':pd.__version__, 'scipy':scipy.__version__, 'scikit_learn':sklearn.__version__,
        'synthetic_demo':not bool(args.observations), 'simulation_steps':args.simulate,
        'initial_fit_warnings':first_info['warnings']})
    (destination/'run_settings.json').write_text(json.dumps(info, ensure_ascii=False, indent=2), encoding='utf-8')
    save_plots(initial_observations, first_ranked, history if args.simulate else [], destination)
    print(ranked.head(3).to_string(index=False))
    if info['warnings']:
        print('Model fitting warnings: inspect run_settings.json before using the proposal.')
    print('Saved to', destination.resolve())

if __name__ == '__main__':
    main()