あるケミストの独り言(winchemwinの日記)

ケミスト(化学者)の視点で、面白そうな情報(シミュレーション関係など)を発信

機械学習による予測モデル作成(我が家の電気使用量 その1-1)

 様々なデータの解析に機械学習の手法が取り入れられてきています。python機械学習のライブラリが比較的充実していますので、今回から機械学習によるいくつかのデータの予測モデル作成と検証についての内容を書いてゆきたいと思っています。

 まずは、我が家の電気使用量について行ってみました。我が家は関西電力さんと契約をしているのですが、関西電力さんの「はぴeみる電」のホームページから各自日々の電気使用量をダウンロートできます。そのデータのうち電気使用量を目的変数として、説明変数との関係について機械学習のモデルの作成を検討しました。

 説明変数ですが、なかなか関係が深そうなデータで系統的に集積されているものはないのですが、今回は気象庁で公開されている日々の気象情報との関係について、機械学習で整理できないか検討してみました。気温や湿度はエアコンなどの使用量と関係するので、それなりの相関関係があるのではないかと期待した訳です。

 

 以下はそのコードですが、まずはデータの読み込みと前処理です。それぞれのデータですが、電気使用量は、1ヶ月単位で「はぴeみる電」のホームページからダウンロードできます。

データの読み込みと前処理

 

データは毎日、毎時間の使用量と簡単な気象情報からなっていますが、今回は毎日の電気使用量を目的変数として設定することにしました。データは1ヶ月ごとのデータなので、年間のデータとして予め1年間分のデータをダウンロードして、

集約したファイルを「Year2020_EV.csv」としてエクセルで編集しています。Pandas のデータフレームとしてデータを読みこんでいますが、カラム名が割当られていないのでそれぞれのカラム名を df.renameでわかりやすい名前に割り当てています。

 元のデータの日付ですが、”4/1”のように年のデータが含まれ

ていないので’2020/’を加えた上でpd.to_datatime で日時データに変換しています。また電気使用量はウィークデイと週末とで大きく異なるので、曜日による分類ができるようdt.strftimeで曜日情報を取得した上で、新たにDay/Weekとしてデータフレームに列を追加しています。以上で電気使用量データの前処理は終了です。

 続いて気象データです。気象データも各地の気温、風速、気圧、湿度などがダウンロードできますが、一度にすべてのデータをダウンロードできるわけではないのでそれぞれの個別にダウンロードしたものをエクセルで集約、編集しています。気象データも毎日、毎時間のデータがありますが、今回は毎日のデータを説明変数として利用しています。

 加工した気象データ(Weather2020.csv)はPandas のデータフレームとして読み込んでいます。ここで、電気使用量データと気象データの日

付の並びが異なるので、気象データの方の順番をdf.iloc[::-1]で逆に並び替えて、電気使用量データをあわせる作業を行っています。続いて電気使用量と関係しそうな必要データを取り出して、気象データの前処理作業は終了です。

 最後に電気使用量データと気象データを統合して、機械学習用データ(dfsum)としています。また電気使用量ですが、ウィークデイと週末では大きく異なるので、それぞれ分けて機械学習モデルの作成を行うため、データを分割しています。

  

 以上、説明が長くなりましたが、電気使用量のデータ、気象データの読み込みと前処理までの説明でした。

 次回以降、sklearn による機械学習の実施と検証について紹介したいと思います。

Psi4による量子化学計算ー構造最適化(座標の確認)

前回までの記事で水分子の構造最適化の実施について、紹介しました。

通常は最適化計算後、座標データや構造の確認を行うかと思います。

以下は前回紹介したコードですが、計算のアウトプットファイルは"h2o-opt.txt"に保存されていて、このファイルを確認することで計算の詳細、過程や座標データを確認することができます。

水分子の構造最適化ー(前回のもの)

 以下はoutput ファイル内に記載された内容の一部になります。最適化

された座標データは一番したの欄で確認できます。ただ、いちいち

アウトプットファイルを確認するのも手間がかかりますので、座標データ

を確認できるコードを加えてみました。

 

水分子の構造最適化ーアウトプットファイルの内容

 以下は前回のコードにxyz座標データを取り出すコード(save_string_xyz_file( ) ) を加えたものになります。座標データはxyz形式としてファイル保存するとともに、print 文で表示もしています。座標データのみのファイルがあると別の計算を実施したりする際にも便利ですよね。

 

水分子の構造最適化ー座標データの表示

 上記のコードをjupyter notebook上で実施すると以下のような結果が表示され、最適化後の座標が確認できます。座標だけではどのように最適化されているかはイメージしにくいかもしれませんが、小さな分子であればなんとなくそれぞれの原子がどのように位置しているかは想像できるかもしれませんね。

 

表示された座標データ

 次回は、最適化構造の描画について紹介したいと思います。

Psi4による量子化学計算ーインプットファイルの作成 その2

 

 前回までの記事で、座標データの取得まで説明してきましたので、ここからは、実際の計算用コードの作成に移ってゆきたいと思います。以下は水分子の構造最適化のコードの例です。

f:id:winchemwin:20220107161852j:plain

H2Oの構造最適化の参考コード1


 最初にpsi4のライブラリーをインポート(import psi4)した後、計算に使用するスレッド数とメモリーを設定します(psi4.set_num_threads(nthread=##), psi4.set_memory('###MB'))。私の場合は低スペックのPCで行っていたので上記のコードに記載のスレッド、メモリーですが、PCのスペックに合わせて適宜設定して頂いたら良いかと思います。

 次に化合物の構造情報を記載していきます。座標は前回の記事で記載したxyz形式の座標データを使っています(h2o=psi4.geometry(""" 座標データ"""))。

 続いて、計算結果を保存するためのファイルとしてアウトプットファイルの設定も行った上で(psi4.set_output_file('ファイル名’)、計算手法(今回は’hf/3-21g')、対象分子を設定した上で構造最適化を行うコード(psi4.optimize(計算手法、対象分子)となっています。

f:id:winchemwin:20220107164151j:plain

H2Oの構造最適化結果その1

 

 計算を実行すると最初に設定されたスレッド数とメモリーが表示されます。しばらくして構造最適化に成功するとOptimization completeの表示が現れるとともに、安定化された構造のエネルギーが表示されて計算が終了します。別の化合物を計算する場合は化合物名、座標データ、ファイル名、計算手法などを変えてコードを作成し、計算させることで同様に構造最適化計算を実施できます。

 以上、シンプルに構造最適化する手法を紹介させていただきましたが、次回以降

は最適化構造の座標の表示や化学構造の描画について、紹介したいと思います。

 

 

pythonの活用例について(fit ファイルの利用 その15:Cateye社 Avventuraのデータからーアプリの作成)

 前回までの記事で、fit ファイルをもとにpython で様々な解析を行ってきた例を紹介してきました。それぞれの解析は紹介したコードを参考に処理してもらえればできるかと思いますが、やはりいちいちコードを書いて実行するのは・・・ということでtkinter を使って簡単なアプリを作ってみました。下記が作成したアプリを立ち上げた画面になります。

 使い方ですが、「解析用のファイルを選択」と書かれたところの下にある「選択」ボタンを押して、解析用のfit ファイルを選択します。選択したファイルが上の空欄のところにディレクトリも含め表示されたらOKです。

 次に選択したfitファイルの解析範囲を設定するため(解析範囲を距離で設定)、その下の「解析範囲を設定」ボタンをおします。

アプリ メイン画面

 次に選択したfitファイルの解析範囲を設定するため(解析範囲を距離で設定)、その下の「解析範囲を設定」ボタンをおします。すると別画面に標高/距離のグラフが表示されます。(下図)グラフを参考にしながら、「解析開始点」、「解析終了点」を入力して解析したい範囲を指定します。数値を入力し終えたら、右下の「実行」ボタンを押して、解析範囲の設定を実行します。

解析範囲の設定画面

 解析範囲の設定が終わると、下図のように解析範囲がメイン画面に表示されます。

 

設定範囲の表示(解析範囲設定の結果)

 次にメイン画面の「作成するグラフ、データ等を選択してください」書かれている下のプルダウンメニューから行いたい処理を選択します。選択した処理はその下の「実行」ボタンを押すことで実行されます。解析結果は別画面で例えば下図(いくつかあります)のように表示されます。解析結果の表示画面で表示画面の拡大縮小や画面の保存も行うこともできます。

距離/速度、勾配のグラフ

 

ヒストグラム(勾配)

 

バイオリンプロット(勾配)

激坂ファクター

 読み込んだファイルはメイン画面の一番下の「データ保存」のボタンを押すことでcsv形式のファイルに保存することも可能です。

 シンプルなアプリですが、よかったら試してみてください。

 

 尚、このアプリのコード(英語版)は下記のgithub にMITライセンスで公開していますので、そちらから利用してみてもらったらと思います。

 

github.com

 

Psi4による量子化学計算ーインプットファイルの作成 1

デモのプログラムでPsi4の動作が確認できたら、いよいよ実際の掲載に移っていきます。

Psi4では化合物の座標情報として、デカルト座標, Z-matrix 系 のインプットの入力が可能です。ここではデカルト座標系での入力で計算を実施してゆきたいと思います。化合物の構造から座標データを得るには、モデリングソフトなどを利用するのが簡単です。市販の有料ソフトウェアを使って作成することは可能ですが、無料で簡単に行うとなるとWeb上で構造式描画を行えるMolviewがひとつの選択肢になります。

 Molvewについては下記のリンクからアクセスすることができます。

 

 MolView

 

 リンク先にアクセスすると下記の初期画面が現れます。利用規約などを確認した上で”Close"をクリックすると次の初期入力画面(カフェインの構造がデフォルトで入力されている)が現れます。

f:id:winchemwin:20211212154453j:plain

Mol View初期画面

 

 

 

f:id:winchemwin:20211212155543j:plain

初期入力画面

 MolViewの使い方ですが、初期入力画面からカフェインの構造を消去した後に、作成したい構造を書いていきます。下記は入力画面ですが、左側に必要な構造のパーツがあり、右側に主要元素が並んでいます。化学系の描画ソフトを使ったことのある方なら、それほど違和感なく書きたい構造をかけるのではないかと思います。構造の書き方はこのブログでは詳しくは説明はしませんが、それほど難しいソフトではないので触ってなれてみてもらえたらと思います。構造がかけたら、3Dデータを作成するため、上部のメニューバーにある2D to 3Dをクリックします(これが大事なポイントです)。すると右側の画面に描画した分子の3Dの分子モデルが表示されます。

f:id:winchemwin:20211218040435j:plain

入力画面

  下記は「酢酸エチル」を入力して、3D表記させた例になります。

f:id:winchemwin:20211218042559j:plain

入力例(酢酸エチル)

 構造を作成できたら、Mol ファイルへの保存を行います。ToolのプルダウンメニューからExport - MOL file を選択します。

 

f:id:winchemwin:20211218043037j:plain

Mol file での保存 1

 保存を行うとソフトウェアの下部に保存したファイル名が表示されます。表示されたファイルをダウンロードして自分のPCにファイルを保存することで、作成した構造のMol fileを取得することができます。

 

f:id:winchemwin:20211218043628j:plain

Mol file での保存 2

 以上でMol file を作成することができましたので、次にインプットファイルの形式に対応したデカルト座標のデータを得るため別のソフトウェアで処理を行います。単純にファイル形式の変更ですので、mol ファイルからxyz形式に変更できる手法であればここで紹介する手法以外(openbabelなどの利用)でも問題はないかと思います。ここでは無料利用可能なAvogadroを使った方法を紹介します。Avogadroについては下記リンクからアクセスすることができます。

 

Avogadro - Free cross-platform molecular editor - Avogadro

 

内容確認の上、使用OSにあったプログラムをダウンロードしてもらえればと思います。

f:id:winchemwin:20211226035936j:plain

Avogadro ホームページ画面

導入後、ソフトウェアを立ち上げると下記のような画面が現れます。

 

f:id:winchemwin:20211226062805j:plain

Avogadro 立ち上げ画面

立ち上げ後、作成したmolファイルを読み込み、改めてxyzとして保存します。

f:id:winchemwin:20211226063314j:plain

保存画面

保存したファイルを確認すると下記のようにデカルト座標系のデータとなっていますので、このデータを用いてPsi4のインプットファイルを作成してゆきます。

f:id:winchemwin:20211226063711j:plain

座標データ(H2Oの場合) 

実際のPsi4のインプットの作成は次の記事で紹介します。

pythonの活用例について(fit ファイルの利用 その14:Cateye社 Avventuraのデータからーヒルクライムデータ-坂データの比較)

これまでの記事でfitファイルのデータから、様々な解析を行って来ました。特にヒルクライムに関連して、勾配分布や関連グラフ、自作の激坂ファクターの例などを紹介してきました。今回はそれらの集大成(すこし大げさですね)として、いくつかの実際の坂データの比較を行ってみたいと思います。

 以下の2つのデータは著者の近隣の坂のデータを解析した結果になります。一般的によく使われる平均勾配と距離で比較すると「A」と「B]は似たような数値となっていますが、実際に実走すると、明らかに「A」の坂の方がきつく感じます。

 

f:id:winchemwin:20220320105912p:plain

坂データの比較

 fit ファイルの解析結果で勾配分布を見ると「A」の方が10%前後の勾配の坂の割合が多く、「B」と比べると坂全体としてきつい勾配の割合が高いことがわかります。またバリオプロットでも同じような傾向が確認できるのがわかるかと思います。激坂ファクターでもきついの勾配の割合を強調したファクター4、5で「A」の方が大きくなっています。

 このように坂の視覚化という意味ではfitファイルの解析結果は使えるデータもあるかと思いますので、今回までの記事が皆さんの参考になれば幸いです。

 

pythonの活用例について(fit ファイルの利用 その13:Cateye社 Avventuraのデータからーヒルクライムデータ-激坂ファクター)

 前回までにヒルクライムに参考になる図として、勾配に関する箱ひげ図やバイオリンプロットについて紹介してきました。これらの図は坂全体の様子を把握するには非常にイメージしやすいのですが、実際にいろんな坂を比較して「どのくらいきついの」ということがすぐに分かる数字があればわかりやすいかと思いますが、いかがでしょうか。一般的には平均勾配や最大斜度の数字から判断するかと思うのですが、激坂度を一言で示す「ファクター」があればよりわかりやすいのでは?と思い、「激坂ファクター」なるものを作ってみました。

 以下はそのコードですが、基本的に勾配ごと(3-6%、6.5-9%、9.5-12%、12.5-15%、15.5-18%)の距離(Sumdist)を算出した上で、勾配大きい距離がどのくらいあるかでファクターを算出しています。

f:id:winchemwin:20220312042714j:plain

激坂ファクター コード その1(距離の算出)

f:id:winchemwin:20220312042738j:plain

激坂ファクター コード その2 (距離の算出)

 ファクターですが、いくつか算出しています。ファクター1は全体の距離に対する「15.5-18%」の割合、ファクター2は「12.5-15%」の割合、ファクター3は「9.5-12%」の割合を示しています。ファクター1-3は単純にそれぞれの勾配の比率になっていますので、それとは別にファクター4、5は激坂を強調するために勾配の大きい部分に大きな係数をかけて計算してみました。

 具体的にはファクター4ではそれぞれ「9.5-12%」✗1、「12.5-15%」✗2、「15.5-18%」✗3した距離に対する全体の距離、ファクター5ではそれぞれ「3-6%」✗0.5、「6.5-9%」✗1、「9.5-12%」✗1.5、「12.5-15%」✗3、「15.5-18%」✗5した距離に対する全体の距離で算出しています。

f:id:winchemwin:20220312042804j:plain

激坂ファクター コード その3(ファクターの計算)

 以下は、ある坂について激坂ファクターを算出してみた結果になります。数字だけみてもピンとこないかもしれませんが、いろんな坂で比較してみると「激坂ファクター5」が「1」を超えると個人的には「きつーい」坂に分類できているように思います。興味があれば皆さんのデータでも計算されてみてもいかがでしょうか。

f:id:winchemwin:20220312043302j:plain

激坂ファクター?? 算出結果

  次回はよりイメージしやすいようにいくつかの坂でバイオリンプロットと激坂ファクターを計算してみた結果について紹介したいと思います。