pymatgenでVASPの入力ファイルを作る

pymatgen

pymatgenを使うとVASPの入力ファイルを簡単に作ることができるが、公式ドキュメントにはまとまった記述がなかったので自分用に備忘録を残しておく。 正直、pymatgenのAPIドキュメントを読めばこの記事を読む必要はないと思う。

POTCAR

http://pymatgen.org/installation.html#potcar-setup に従ってpymatgenが読むPOTCARの配置場所を設定する。

まず、例えばVASPに付属しているPOTCARのディレクトリ(<EXTRACTED_VASP_POTCAR>)を<MY_PSP>以下にpymatgenが読み込める形式でコピーする:

pmg config -p <EXTRACTED_VASP_POTCAR> <MY_PSP>

次に<MY_PSP>をpymatgenが読めるようにする:

pmg config --add PMG_VASP_PSP_DIR <MY_PSP>

実は上のコマンドを実行することによって、ホームディレクトリ直下に以下の内容の.pmgrc.yamlを作成されている。

以上の設定でpymatgenからPOTCARを読み込めるようになる。 POTCARはpymatgen.io.vasp.inputs.Potcarで作ることができる。

POSCAR

pymatgen.core.Structureからpymatgen.io.vasp.inputs.Poscarで作ることができる。 既存のPOSCARからpymatgenのPoscarオブジェクトを作りたい場合はpymatgen.io.vasp.inputs.Poscar.from_fileを使えば良い。

一般にPOTCARにおける原子種の並びとPOSCARにおける原子種の並びは同じになっていなければならない(http://cms.mpi.univie.ac.at/wiki/index.php/POTCAR)。 そのため、pymatgen.io.vasp.inputs.Potcarに渡す原子種(symbols)とpymatgen.core.Structureの原子種(species)は同じになっていなければならない。 具体的には

potcar = pymatgen.io.vasp.inputs.Potcar(symbols)
structure = pymatgen.core.Structure(lattice, species, coords)
poscar = pymatgen.io.vasp.inputs.Poscar(structure)

みたいなときに

symbols == [a[0] for a in itertools.groupby(species)]

が成り立っていなければならない。

INCAR

INCARタグをキーとするdictからpymatgen.io.vasp.inputs.Incarで作ることができる。 既存のINCARからpymatgenのIncarオブジェクトを作りたい場合はpymatgen.io.vasp.inputs.Incar.from_fileを使えば良い。

KPOINTS

APIドキュメントを読めば分かるので省略。

例

VASPの公式チュートリアルのfcc Si を計算する場合の例を示す。 交換相関エネルギー汎関数にはPBEを用いた。

import numpy as np

from pymatgen.core import Structure, Lattice
from pymatgen.io.vasp.inputs import Incar, Kpoints, Poscar, Potcar, VaspInput


def fccSi(a):
    cell = a * np.array([[0., 0.5, 0.5],
                         [0.5, 0., 0.5],
                         [0.5, 0.5, 0.]])
    coords = np.array([[0., 0., 0.]])

    struct = Structure(Lattice(cell), ['Si'], coords)
    return struct


def write_input(a, output_dir='.'):
    # poscar
    structure = fccSi(a)
    comment = None
    poscar = Poscar(structure, comment)

    # incar
    params = {
        'System': 'fcc Si',
        # not read WAVECAR, and start from superposition of atomic charge density
        'ISTART': 0, 'ICHARG': 2,
        'ENCUT': 240,
        'ISMEAR': 0, 'SIGMA': 0.1
    }
    incar = Incar(params)

    # kpoints
    kpoints = Kpoints(num_kpts=0, style="Monkhorst", kpts=((11, 11, 11), ),
                      kpts_shift=(0, 0, 0))

    # potcar
    potcar = Potcar(['Si'], 'PBE')

    # write input files
    vaspinput = VaspInput(incar, kpoints, poscar, potcar)
    vaspinput.write_input(output_dir)


if __name__ == '__main__':
    for a in np.linspace(3.5, 4.3, num=9):
        output_dir = "fccSi_{}".format(int(10 * a))
        write_input(a, output_dir)

周期カーネルに対するFeature Map

ガウス過程回帰の課題のひとつとしてスケーラビリティがあります。
データ数をNとしてfitting の時間計算量 O(N^3)、predict の時間計算量が O(N^2)であることから大量のデータを扱いたいときは何らかの近似が必要です。

そのような近似手法のひとつとして、カーネルを有限次元の特徴ベクトルの内積で近似してベイズ線形回帰に帰着するものがあります。
特徴量ベクトルの次元をDとして、ベイズ線形回帰のfitting の時間計算量は O(D^3) 、predict の時間計算量は O(D)なので特徴量ベクトルの次元がデータ数よりも小さければより高速に計算ができることになります。

与えられたカーネルに対応する特徴量ベクトルを作る方法としてRandom Fourier Features
があります。*1
Gaussian カーネルやLaplacian カーネルについてはRandom Feature Map で特徴量ベクトルを得られますが、周期カーネルにこの方法を適用するのは難しそうでした。

そこでいろいろ調べたところ、周期カーネルを近似する論文*2を見つけたので概要と実装を備忘録としてまとめておきます。

(追記: 一部の数式が上手く表示されないようなので画像で置き換えました。何もわからない。)

Random Fourier Features

そもそもRandom Fourier Features がどのようなアルゴリズムなのか簡単に述べておきます。
\mathbb{R}^{D}上のstationaryなカーネル k(\mathbf{x}, \mathbf{y}) = k(\mathbf{x} - \mathbf{y})がpositive definite であることと、ある有限非負測度のフーリエ変換で書けることが同値である(Bochner の定理)ことを使います。

Bocher の定理より、適当に定数倍したpositive definite なカーネルはある確率密度関数 p(\mathbf{\omega})のフーリエ変換で書けます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&= \int_{\mathbb{R}^{D}} p(\mathbf{\omega}) e^{i (\mathbf{x} - \mathbf{y})^{T} \mathbf{\omega}} d\mathbf{\omega} \\
&= \mathbb{E}_{\mathbf{\omega}, b} \left[ z_{\mathbf{\omega}}(\mathbf{x}) z_{\mathbf{\omega}}(\mathbf{y})\right]
\quad (\mbox{where} \quad z_{\mathbf{\omega}}(\mathbf{x}) := \sqrt{2} \cos(\mathbf{x}^T \mathbf{\omega} + b))
\end{align}
最後の式において、 \mathbf{\omega}は p(\mathbf{\omega})に従い、bは [0, 2\pi]の一様分布に従います。

\mathbf{\omega}とbをそれぞれD個ずつサンプリングして
\begin{align}
\mathbf{z}(\mathbf{x}) := \sqrt{\frac{2}{D}} \left( \cos(\mathbf{x}^T \mathbf{\omega}_{1} + b_{1}), \cdots, \cos(\mathbf{x}^T \mathbf{\omega}_{D} + b_{D})\right)^{T}
\end{align}
とすると、カーネルをモンテカルロ積分の要領で近似できます。
\begin{align}
k(\mathbf{x} - \mathbf{y})
&\simeq \frac{1}{D} \sum_{j = 1}^{D} z_{\mathbf{\omega}_{j}}(\mathbf{x}) z_{\mathbf{\omega}_{j}}(\mathbf{y}) \\
&= \mathbf{z}(\mathbf{x})^{T} \mathbf{z}(\mathbf{y})
\end{align}

例えばGaussian カーネルは正規分布をフーリエ変換で書けるのでこの方法で対応する特徴量ベクトルを得ることができます。
これは既にsklearn に実装されています。*3

Fourier Features for Periodic Kernels

次に周期カーネルを考えます。
具体的には次のようなカーネルです。*4
\begin{align}
k_{perSE}(t, t') = k_{perSE}(t, t')
&= \exp \left( -2 \left( \frac{\sin (\omega_{0} (t - t'))}{l} \right)^{2} \right) \\
&= \exp \left( \frac{\cos \omega_{0} (t - t') - 1}{l^{2}} \right)
\end{align}
ここで\omega_{0} = \frac{2\pi}{T}で、このカーネルは周期Tの関数を説明したい時などに使われます。

Random Fourier Features で特徴量ベクトルを作るなら、周期カーネルを逆フーリエ変換して得られる確率密度関数に従って周波数をサンプリングすることになりますが、どうみても簡単にサンプリングができる確率分布を得られそうにありません。
そこで、周期カーネルをフーリエ級数展開して近似することを考えます。

周期カーネルのフーリエ級数展開を次のように書きます。
f:id:lan496:20180716013606p:plain

k_{perSE}(\tau)は偶関数であることと、第一種変形ベッセル関数の積分表示I_{n}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x \cos \theta} \cos n \theta\, d\theta \, (n \in \mathbb{Z})
*5を使うと、 c_{k}は次のように計算できます。
f:id:lan496:20180716013619p:plain

フーリエ級数の第K項以降を切り捨てて近似します。
f:id:lan496:20180716013630p:plain
よって、周期カーネルを2K次元の特徴量ベクトル \mathbf{\phi}(t)で近似できることがわかりました。

実装とテスト

論文の公式実装が見当たらなかったので自分で書きました。*6

論文ではランダムに2000個の入力に対して、周期T=2の周期カーネルでグラム行列 \mathbf{K}を計算したものと、第C項までを使って近似で得られる行列 \tilde{\mathbf{K}}のフロベニウスノルムでの誤差 \frac{|| \tilde{\mathbf{K}} - \mathbf{K}||_{F}}{|| \mathbf{K} ||_{F}} を測っていたので、自分の実装で検証してみました。

特徴量ベクトルの次元とハイパーパラメターlength_scale(l)を変えた時の誤差は下図のとおりです。
図の左上部分の誤差が元論文と異なりますが、傾向としては同じ結果が得られました。
f:id:lan496:20180716011413p:plain

下図はlength_scale(l)=0.2の場合にカーネルの値をプロットしたものです。
フーリエ級数展開の16次まで取り込めば、元の周期カーネルをほぼ再現できていることがわかります。
f:id:lan496:20180716012425p:plain

*1:https://people.eecs.berkeley.edu/~brecht/papers/07.rah.rec.nips.pdf

*2:http://web.it.usyd.edu.au/~framos/Publications_files/Tompkins-Ramos.pdf

*3:sklearn.kernel_approximation.RBFSampler — scikit-learn 0.19.1 documentation

*4:sklearn.gaussian_process.kernels.ExpSineSquared — scikit-learn 0.19.1 documentation

*5:Modified Bessel Function of the First Kind -- from Wolfram MathWorld

*6:feature_map.py · GitHub

化合物の構造を表現する記述子

無機化合物の化学式とユニットセルの情報及び原子の座標から物性値(結合エネルギー・バンドギャップ)を予測するKaggleコンペに参加していた。

Nomad2018 Predicting Transparent Conductors | Kaggle

feature engineering の過程で分子および無機結晶の構造から特徴量抽出する手法に関する論文をいろいろ読んだので、以下にまとめておく。 基本的に記述子の作り方の部分しか読んでいないので、詳しくは論文を参照。

分子に対する手法

Coulomb Matrix (CM)

ここで Z_{I}は原子Iの形式電荷、 r_{IJ}は原子Iと原子Jの距離。

  • 分子を構成する原子Iと原子Jの間の反発力を表している。 対角成分については、原子の自己エネルギーを多項式でフィッティングしたもの。

  • Coulomb Matrix M の固有値を絶対値の大きい順で並べ、適当に0でpadding して次元を揃えたものを記述子として使う(分子ごとに原子数が異なるので記述子の次元を揃える必要がある)。

  • 分子の回転・並進移動の元で記述子は不変。

Bag of Bonds (BoB)

Bonds, angles, machine learning (BAML)

Extended Connectivity Fingerprints (ECFP)

Molecular atomic radial angular distribution (MARAD)

Histogram of distances, angles and dihedral angles (HDAD)

Fourier series of radial distance distributions

Molecular Graph (MG)

無機結晶に対する手法

分子に対する手法ほど進展していない印象。 分子の場合は並進移動と回転に対する対称性を考えればよかったが、結晶の場合は更にユニットセルの選び方に関して不変であるべきなのでより難しそう。 crystal net *2を一意に表す方法があればいいのだが、そもそもグラフ同型判定問題はNPに属しているという厳しさがある。

Partial radial distribution function (PRDF)

Atom centered symmetry functions

f_{c}(R) = \begin{cases} \frac{1}{2} \left( 1 + \cos \left( \frac{\pi R}{R_{cutoff}} \right) \right) \quad (R &lt; R_{cutoff}) \\ 0 \quad (R > R_{cutoff}) \end{cases}

  • 各原子 iに対して次のような動径分布関数を考える。 \eta, R_{s}をいろいろ変えて記述子とする。

G_{i}^{rad} (\eta, R_{s})= \sum_{j} e^{-\eta (R_{ij} - R_{s})^{2} } f_{c} (R_{ij})

  • 各原子 iに対して次のような角度分布関数を考える。 \zeta, \etaをいろいろ変えて記述子とする。

G_{i}^{ang} (\zeta, \eta) = \frac{1}{2^{\zeta - 1}} \sum_{j, k \neq i} \left( 1 \pm \cos \theta_{ijk} \right)^{\zeta} \, e^{-\eta (R_{ij}^{2} + R_{jk}^{2} + R_{ki}^{2})} f_{c} (R_{ij}) f_{c} (R_{jk}) f_{c} (R_{ki})

Bond orientational order parameter (BOP, BOOP)

Bispectrum

F-Fingerprint

F_{AB}(R) = \sum_{A_{i} \in unit cell} \sum_{B_{j} s.t. R_{ij} \lt R_{max}} \frac{\delta (R - R_{ij}) }{4 \pi R_{ij}^{2} \frac{N_{A} N_{B}}{V} \Delta} - 1

  • ここで、 R_{ij}は原子 A_{i}と原子 B_{j}の距離、 N_{A}, N_{B}はユニットセル内の原子A,Bの個数、 Vはユニットセルの体積、 \Deltaはこの関数をヒストグラムにするビン幅、 \delta (\cdot )はgaussian-smeared delta function

Ewald Sum Matrix

  • [1503.07406] Crystal Structure Representations for Machine Learning Models of Formation Energies
  • Coulomb Matrix における M_{IJ} をサイトIを占めるすべての原子とサイトJを占めるすべての原子との間の静電ポテンシャルの和に修正。
  • Ewald法を用いて静電ポテンシャルの総和を計算する。
  • pymatgen にEwald法を用いて結晶の静電ポテンシャルの総和を計算する関数があるが、計算する式が少し違う。ただEwald Sum Matrix を実際に実装するときに大いに参考にはなる。

Extended Coulomb-like Matrix

Sine Matrix

Crystal Graph (CG)

*1: 二面角 - Wikipedia

*2: Periodic graph (crystallography) - Wikipedia

マルチカノニカル法で個数推定

はじめに

去年の11,12月にマルコフ連鎖モンテカルロ法(MCMC)の集中講義を受けてきました。 講義の中で紹介されたマルチカノニカル法を応用してレアイベントをサンプリングする手法が特に面白かったので、復習も兼ねて実装してみました。

マルチカノニカル法

まずレアイベントサンプリング用のマルチカノニカル法を説明します。

状態xの確率分布が\pi(x)が与えられているときに、統計量\xi (x)の分布が知りたいとします。 単純にxをサンプリングしてもいいのですが、統計量の周辺分布P(\xi)に極めて確率の小さい区間があるとその部分をうまくサンプリングすることができません。

この問題を回避するために適当な重みQ(\xi )を用意して、\pi (x)の代わりに重みQ(\xi (x)) \pi(x) に従ってxをサンプリングすることを考えます。 これは普通のメトロポリス・ヘイスティングス法で実現できます。 すなわち、状態xから状態x' に遷移する確率を

\mathrm{min}(1, \frac{Q(\xi (x')) \pi(x') p(x \to x')}{Q(\xi (x)) \pi(x) p(x' \to x)})

とすればいいわけです。

重みQ(\xi )として、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)の逆数を用いるのがマルチカノニカル法です。

重み \tilde{P}(\xi (x))^{-1} \pi (x)に従って xをサンプリングしたのときの統計量 \xi の周辺分布 P'(\xi )は

P'(\xi ) \propto P(\xi) \tilde{P}(\xi )^{-1}

となって、この分布はほぼ平坦であることがわかります。 平坦な周辺分布 P'(\xi )からサンプリングすることで、普通にサンプリングすると確率0と推定されてしまうようなレアイベントにもMCMCで訪れられるのがポイントです。

Wang-Landau 法

残る問題は、求めたい周辺分布P(\xi)に近い分布\tilde{P}(\xi)をどのように得るかです。

周辺分布P(\xi)を大まかに求める方法としてWang-Landau 法というものがあります。 詳細はwikiの擬似コードを読むのが一番分かりやすいと思います(https://en.wikipedia.org/wiki/Wang_and_Landau_algorithm)。

大まかに言うと、 \xiのヒストグラムが平坦になるように出現回数の低いところの重みを逐次大きくしていき、ある程度平坦になったらより細かく重みを調整していく感じです。

演習問題

具体例として、次のような数え上げ問題を考えます。

L×Lの正方格子において、角と向かい側の角とをつなぐ自己交差のない経路の個数を求めよ。

文章だと分かりづらいですが、要は↓の動画で扱われている問題です。 www.youtube.com

一辺の長さがLのとき、辺は全部で2L(L+1)本あるので、各辺について使うか否かを考えると辺の配置は全部で2^{2L(L+1)}通りです。探索する配置の数が指数関数で増大するので、愚直に全通りを試していてはすぐに計算が終わらなくなってしまいます。

ただ、この問題には効率的なアルゴリズムがあり、L=26までの厳密解が求まっているようです(http://oeis.org/A007764?language=japanese)。 なので個数推定をする意義はあまりないのですが、演習問題を解いているということでそこら辺は気にせずいきます。

この問題をマルチカノニカル法で扱うために、正方格子の各辺について使うか否か決めたものC (以下これを「状態C 」と呼ぶ)に対してエネルギーを次式で定めます。

\displaystyle { E(C) = \sum_{i : \mbox{左上と右下以外の頂点}} |\mbox{(iの次数)} - 2| + \sum_{i : \mbox{左上と右下の頂点}} |\mbox{(iの次数)} - 1| + (\mbox{2個以上の頂点からなる連結要素の数)} }

ここで、E(C) = 0であることと状態Cが左上と右下の頂点をつなぐ自己交差のない経路であることが同値になっています。 先のマルチカノニカル法を用いてエネルギーの分布P(E)を求めて、P(E=0)に状態の総数2^{2L(L+1)}をかければ今求めたい経路の数を推定することができます。

具体的な実装はこちら。 github.com

結果

以上の方法で自己交差のない経路の数を100回推定し、その平均と不偏分散の平方根を真値と比較したものが下図です。 有効数字2桁くらいの精度が得ることができました。

ちなみにL=6が手元のMacBook Air で4時間くらいかかっています。 L=6の状態の数が 2^{2 \times 6 \times 7} \sim 3 \times 10^{50}で、その中で自己交差のない経路になっているものはおよそ 6.4 \times 10^{42}分の1しかないことを考えると、これほどレアな事象をそれなりの精度でサンプリングできるのはすごいと思えます。

L 推定値 誤差 真値
3 185.072 2.6228 184
4 8241.03 275.432 8512
5 1.2242e+06 75237.9 1262816
6 5.58901e+08 6.03206e+07 575780564

参考文献

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125062 https://sites.google.com/site/iwanamidatascience/vol2/number-place http://www.springer.com/us/book/9783642031625

AOJ 1320 City Merger

AOJ

問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1320)

AOJ2022(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2022)とほとんど同じ問題。他の都市名の部分文字列になっている都市名を予め取り除いてからbitDPをする。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define ALL(x) (x).begin(), (x).end()

typedef long long ll;
typedef long double ld;

const int INF = 1e9;
const ld EPS = 1e-8;

bool is_contained(string s, string t) {
  for(int i = 0; i + s.length() - 1 < t.length(); ++i) {
    if(s == t.substr(i, s.length())) return true;
  }
  return false;
}

int d[16][16];
int memo[16][1<<16];

int dp(int n, int k, int S) {
  //cout << k << " " << S << endl;
  if(memo[k][S] != INF) return memo[k][S];
  if(S == (1<<k)) return d[k][k];
  int res = INF;
  REP(i, n){
    if(i == k) continue;
    if(S & (1<<i)){
      res = min(res, d[k][i] + dp(n, i, S - (1<<k)));
    }
  }
  //cout << k << " " << S << ": " << res << endl;
  return memo[k][S] = res;
}

int main(){
  int n;
  while(cin >> n && n) {
    vector<string> cities(n);
    REP(i,n) cin >> cities[i];
    vector<string> tmp;
    REP(i,n){
      bool flag = true;
      REP(j,n){
        if(i == j) continue;
        if(is_contained(cities[i], cities[j])){
          flag = false;
          break;
        }
      }
      if(flag) tmp.push_back(cities[i]);
    }
    cities = tmp;
    n = cities.size();

    REP(i,16)REP(j,16) d[i][j] = INF;
    REP(i,n)REP(j,n){
      string s = cities[i];
      string t = cities[j];
      bool flag = true;
      for(int k = max(0, (int)(s.length() - t.length())); k < s.length(); ++k){
        if(s.substr(k) == t.substr(0, s.length() - k)){
          d[i][j] = k;
          flag = false;
          break;
        }
      }
      if(flag) d[i][j] = s.length();
    }
    REP(i,n) d[i][i] = cities[i].length();
    
    int res = INF;
    REP(i,16)REP(j,1<<16) memo[i][j] = INF;
    REP(i,n){
      int tmp = dp(n, i, (1<<n) - 1);
      res = min(res, tmp);
    }
    cout << res << endl;
  }
  return 0;
}

AOJ 2534 Dictionary

AOJ

問題文(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2534)

まず、与えられた文字列の組(string_i, string_j)(i < j)が辞書式順序を満たすようにアルファベットの間に順序を入れていく。アルファベットx, yに順序が入っていることは有向グラフに(x, y)の辺を張ることで表現する。文字列の長さが違うときは'@'(='a' - 1)で埋めておいた。

このグラフの構築段階で(i,j)に辺があるのに(j,i)にも辺が張られているば矛盾しているのでfalseを返す。また、'@'は最小であるはずなのでそれが満たされていない時もfalseを返す。

'@'を始点として辺の張られているアルファベットに深さ優先探索で移動していくとき、もしも矛盾がなければ高々深さは26にしかならないから、それ以上の深さを探索しようとしていたらfalseを返す。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)

bool dfs(int d[][27], int k, int cnt) {
  if(cnt > 30) return false;
  bool flag = true;
  REP(j,27) if(d[k][j] == 1) flag = false;
  if(flag) return true;

  bool res = true;
  REP(j,27){
    if(d[k][j] == 1){
      res = res && dfs(d, j, cnt + 1);
    }
  }
  return res;
}

bool solve(vector<string> &s) {
  int n = s.size();
  int d[27][27] = {0};
  REP(i,n)for(int j = i + 1; j < n; ++j) {
    int pos = 0;
    while(pos < s[i].length() && s[i][pos] == s[j][pos]) pos++;
    if(pos == s[i].length()) continue;
    if(d[s[j][pos]-'a'+1][s[i][pos]-'a'+1] == 1) return false;
    d[s[i][pos]-'a'+1][s[j][pos]-'a'+1] = 1;
  }
  REP(i,27) if(d[i][0] == 1) return false;
  
  return dfs(d, 0, 0);
}

int main(){
  int n;
  while(cin >> n && n) {
    vector<string> s(n);
    REP(i,n) cin >> s[i];
    int lmx = 0;
    REP(i,n) lmx = max(lmx, (int)s[i].length());
    REP(i,n)REP(j,lmx-s[i].length()) s[i] += 'a' - 1;
    if(solve(s)) cout << "yes" << endl;
    else cout << "no" << endl;
  }
  return 0;
}

Codeforces Round #341 (Div. 2) D. Rat Kwesh and Cheese

Codeforces

問題文はこちら(Problem - D - Codeforces)

まず\( (x^y)^z = (x^z)^y \)であることから\(a_4, a_8, a_{12}\)が候補から外れる。
単純に与えられた式を計算するとオーバーフローしてしまうため、対数をとって比較していく。ただし、\(x < 1\)のとき\( \mathrm{log} x < 0 \)より\( \mathrm{log} \mathrm{log} x\)は(実数上で)存在しない。よってx, y, zの中に1より小さいものが含まれているかどうかで場合分けをする。

(i) x,y,zが全て1より大きい時

\( \{\mathrm{log} \mathrm{log} a_i\}\)の中で最大のものを探す。ただし、最大値を取るものが複数あるときはindexの最も小さいものを選ぶ。

(ii) x, y, zのいずれかが1よりも大きい時

1よりも小さい数が底ならば、今の場合指数は必ず正であるから最終的な値は1よりも小さくなる。逆に1以上の数が底ならば、最終的な値は1以上になる。従って、底が1よりも大きいもののみを(i)と同じように比較すればよい(底が1より小さいものは適当に-INFとかで埋めておく)。

(iii) x,y,zが全て1以下の時

オーバーフローする心配が無いため、\( \mathrm{log} a_i\) の中で最大のものを探せばよい。


コンテストでは\( \mathrm{log} \mathrm{log} x\)が死にうることに気づかず提出してしまった。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define ALL(x) (x).begin(), (x).end()

typedef long long ll;
typedef long double ld;

const int INF = 1e9;
const ld EPS = 1e-8;

bool comp(pair<ld,int> p1, pair<ld,int> p2) {
  if(p1.first != p2.first) return p1.first < p2.first;
  return p1.second > p2.second;
}

int main(){
  double x, y, z;
  cin >> x >> y >> z;
  string ans[] = {"x^y^z", "x^z^y", "(x^y)^z", "y^x^z", "y^z^x", "(y^x)^z", "z^x^y", "z^y^x", "(z^x)^y"};
  vector<pair<ld,int>> val(9);

  if(x > 1.0 || y > 1.0 || z > 1.0) {
    if(x > 1.0) {
      val[0] = make_pair(z * log(y) + log(log(x)), 0);
      val[1] = make_pair(y * log(z) + log(log(x)), 1);
      val[2] = make_pair(log(y * z) + log(log(x)) , 2);
    }else{
      val[0] = make_pair(-INF, 0);
      val[1] = make_pair(-INF, 1);
      val[2] = make_pair(-INF, 2);
    }
    if(y > 1.0) {
      val[3] = make_pair(z * log(x) + log(log(y)), 3);
      val[4] = make_pair(x * log(z) + log(log(y)), 4);
      val[5] = make_pair(log(z * x) + log(log(y)) , 5);
    }else{
      val[3] = make_pair(-INF, 3);
      val[4] = make_pair(-INF, 4);
      val[5] = make_pair(-INF, 5);
    }
    if(z > 1.0) {
      val[6] = make_pair(y * log(x) + log(log(z)), 6);
      val[7] = make_pair(x * log(y) + log(log(z)), 7);
      val[8] = make_pair(log(y * x) + log(log(z)) , 8);
    }else{
      val[6] = make_pair(-INF, 6);
      val[7] = make_pair(-INF, 7);
      val[8] = make_pair(-INF, 8);
    }
  }else{
    val[0] = make_pair(pow(y, z) * log(x), 0);
    val[1] = make_pair(pow(z, y) * log(x), 1);
    val[2] = make_pair(y * z * log(x), 2);
    val[3] = make_pair(pow(x, z) * log(y), 3);
    val[4] = make_pair(pow(z, x) * log(y), 4);
    val[5] = make_pair(x * z * log(y), 5);
    val[6] = make_pair(pow(x, y) * log(z), 6);
    val[7] = make_pair(pow(y, x) * log(z), 7);
    val[8] = make_pair(x * y * log(z), 8);
  }

  sort(ALL(val), comp);
  cout << ans[val[8].second] << endl;
  return 0;
}