pybind11

手は動かしていませんが、忘れそうなのでメモ ( ..)φ

pybind11
https://github.com/pybind/pybind11

pybind11 is a lightweight header-only library that exposes C++ types in Python and vice versa, mainly to create Python bindings of existing C++ code. Its goals and syntax are similar to the excellent Boost.

C++ 側の書き込みは最低限で、私のような初心者でもわかりやすいのがgood!
ネットの書き込みを見ていますと、Python 側の list は C++ 側の vector や list に変換してもらえるとのこと。オーバーヘッドも Cython より小さくなりやすいと評価される方もいらっしゃいました。

本当ならありがたい。まずはこれらから確認してみましょう。。

Infraworks2023でPLATEAU その2

続きです。

1) のテキスト置換を Python で組んだのですが、配布しても使えない方が多い。

ということで、C++で作ってみました。

といっても、C++ は初心者。既存のコードを読んだりコンパイルしたりしたことはありますが、書いたことはありません。久しぶりにプログラミングの入門書を買ってきました。

読んでみると、独学で振り落としてきた重要なパーツがちらほら。半分ほど読んでみましたが、勉強になりました。

で、ネットでコードを拾いながら３日かかって完成（Python なら１時間かからないのに）。特にパスの型が難しかった。意外と、欲しい情報がなかったので、図書を買ったのは正解でした。

https://github.com/T40O0/CityGML-Replace

python も C++ も、バッチファイルを書いて ２) の投影変換と連携させました。ダブルクリックで１）２）を済ませてくれます。できれば C++ で Java まで取り込みたかったのですが、今は実力不足。いずれ整理しましょう。いえ、そのころには PLATAEU で 全都市 FBX が整備されているでしょうね。

Cython

Pandas の公式に速度改善に関するチュートリアルがあります。
https://pandas.pydata.org/pandas-docs/stable/user_guide/enhancingperf.html

CPython ではなく、Cython の利用です。
C で書いたライブラリを Python で読むのではなく、Python で書いた関数を C に渡してコンパイルする形。渡す際に型を定義したり、型チェックを外すことで速くすることが可能になっているようです。これは、使い勝手が良さそう。

プロファイルは動かなかったのですが、それ以外は動きました。遅い処理はこのような簡易な処理方法があることも覚えておきましょう。

C++ & Python

C++ で Pythonのモジュールを作りたい。

初めは Python と C++ の exe をバッチファイルで繋ごうかと思いましたが、ライブラリにすればPython から直接読めます。

チュートリアルがあり、GitHub にテンプレートが公開されていました。それをトレースしましょう。https://docs.microsoft.com/ja-jp/visualstudio/python/working-with-c-cpp-python-in-visual-studio?view=vs-2022

まずは環境設定。
VS2019 は Python3.7 までのサポートのようです。これは古い。VS2019 を削除し、2022に入れなおし。これでデフォが3.9？でしょうか。デフォを3.10にしたかったのですが、入れ替え方がわかりません。
空のPythonプロジェクトを作り、新たに conda 環境を Python=3.10 で作成。これを新しいプロジェクトに対する既定の環境に設定。あらためて空の Python プロジェクトを作り、テンプレートのコードをコピぺしました。が、これではダメ。再度開くと3.9に戻ります。
Google先生に伺ったところ、[全般][インタープリター]で変更できました。
Python のみで動くように初期コードに書き換えた後、3.10で動作確認。OK でした。

次に C++。
既存のプロジェクトとして ’superfastcode’ を追加。VS2017 製でしたので[プロジェクト][ソリューションの再ターゲット]にてアップデート。’superfastcode’ というライブラリとして Python から呼び出せるとのこと。ソースファイル module.cpp は既に追加されている状態です。
プロジェクトのプロパティから、いくつかの設定を変更します。
・[全般][プロジェクトの既定値][構成の種類] dllを指定。
・[C/C++][全般][追加のインクルード ディレクトリ]新たに作成したenv下の include フォルダーを追加。
・[C/C++][コード生成][ランタイムライブラリ]マルチスレッドDLL (/MD
・[リンカー][全般][追加のライブラリディレクトリ] ]新たに作成したenv下の libs フォルダーを追加。

これで一度ビルド。OKです。
’superfastcode2’ も同様に取り込んで設定。採用していた conda 環境に pybind11 をインストールしておきます。念のため、Win11 を再起動してからビルド。
うーん、引っかかります。Python3.6が必要と出ます。どこでしょうか。

結局、原因がわからなかったので’superfastcode2’ はアンロード。
リリース環境でも同じ設定でビルド。で、試算。

Running benchmarks with COUNT = 500000
[tanh(x) for x in d] (Python implementation) took 0.512 seconds
[fast_tanh(x) for x in d] (CPython C++ extension) took 0.123 seconds

やはりCの方が速い。
チュートリアルを見ると pybind11 を使うと ’superfastcode’ より遅くなっていますので、ここまでにしておきましょう。

C++ でライブラリを作り、Python で読み込んで計算するチュートリアルでしたが、分かり難いですね。これ、gcc でも設定できるでしょうか。自信ないですね。

gcc と エンコード

C ++の初心者です。

C、C#とともに何度か入門しましたが、実務ではそれほど用がなかったため、覚える前に辞めてしまう状態が続いています。此度、実務（Windows10）で使うことになり、あらためて入門です。

Visual Studio 2019 を使ってコンパイルしていると、scanf でエラー。

error C4996: 'scanf': This function or variable may be unsafe. Consider using scanf_s instead. To disable deprecation, use _CRT_SECURE_NO_WARNINGS.

そうでした。https://phreeqc.blogspot.com/2021/05/c.html
scanf_s を使うと通りましたが、デバッグでエラー。

test.exe (プロセス *****) は、コード -****** で終了しました。

忘れていますね。scanf_sで文字列を入力する場合には第三引数で最大配列数を指定する必要がありました（引数不足のwarning が出ていましたが、無視していました）。 

引数を指定し、できたexe を実行！
動きました。

表示がすぐに消えてしまうので、「Enterkey を押して終了させたい。」をコピペ。
https://teratail.com/questions/140273
これでOK。

次はgcc
tdm64-gcc-10.3.0-2 を Win11 に入れて、ターミナル起動。gcc -v で動作確認。

>gcc -v
gﾝﾝ spec gpｵﾄ｢ﾜｷB
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/10.3.0/lto-wrapper.exe
・・・
T|[gｳﾄ｢ LTO ｳkASY€: zlib zstd
gcc o[W 10.3.0 (tdm64-1)

文字化けします。
コマンドプロンプト表示をUTF-8に変更

>chcp 65001
Active code page: 65001

>gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/10.3.0/lto-wrapper.exe
・・・
algorithms: zlib zstd
gcc version 10.3.0 (tdm64-1)

OKです。では、コンパイル。コードはUTF-8。

>g++ -o test test.c

できたexeを起動すると、文字化けしています。Win の実行ファイル側では shift-JIS 固定なのでしょう。
実行ファイルでのエンコードをCP932として指定します。

>g++ -o test test.c -fexec-charset=CP932

うまく変換できているように見えます。もう少しテストが必要かな。
ひとまず、環境構築は完了です。

float, double

C#で作成されたソフトと、Pythonで書かれたスクリプト。
同じテキストファイルを読んで計算したのですが、答えが異なりました。

原因は変数の型。

C# は float（32bit） で読み込むように作られていました。Python の float は C の double (64bit)で作られているそうですので、そのまま読み込むと前者の桁が小さくなります。https://docs.python.org/ja/3/library/stdtypes.html

後者を np.float32 で変更すれば良いのでしょうが、そもそも、input ファイルの段階で有効桁数以上の数字を消していなかったのがマズイ。不要な桁を消して読み込めば、一致しました。

手を抜くとダメですね。すぐに帰ってきます。

 

passive 手法

passive 手法のテスト。

*** CC-FJpy ******************************
https://github.com/ColinLii/CC-FJpy
CC-FJpy: A Python Package for seismic ambient noise cross-correlation and the frequency-Bessel transform method
・複数の観測点間の相互相関関数をプロット。
・このコードに限りませんが、基本、長期データが必要です。F-netだと申し込み時に制限（１GB以下、リクエスト数（容量）制限）があるため、難しい。

*** noisepy ******************************
https://noise-python.readthedocs.io/en/latest/
・複数の観測点間の相互相関関数をプロット
・SAC を ASDF に変換して読み込み。
・noise_module が import できないエラー（S0B_to_ASDF.py）。

*** AANT *********************************
http://penglipk.github.io/AANT/
Adaptive ambient noise tomography and its application to the Garlock Fault, southern California
・sac 形式のデータを読み込み、速度を推定し、速度場を逆算してくれます（トモグラフィー）。
・cutsac のソースが含まれておらず中止。
・理想的なツールなのですが、更新が止まっています。非常に惜しい。

*** BIDO *********************************
https://staff.aist.go.jp/ikuo-chou/BIDO/2.0/bidodl.html
・時間領域と周波数領域の結果を比較する目的で実施。
・2.0 では 水平動も利用されていましたが、 3.0 では鉛直動のみ扱うとのこと。
・Ubuntu の binary が入っていましたが、動かなかったのでコンパイルしました。
・15分程度のデータであれば動作。速度が異常。
・1時間程度のデータは読み込み時にエラー。ディメンジョン修正箇所がわからず。

*** FMST *********************************
http://rses.anu.edu.au/~nick/surftomo.html
・トモグラフィーのみのようでした。分散曲線を出した後のステップです。
・速度を事前に用意することで、それらを満たす速度場を逆算。
・マニュアルを読みましたが、今回は手を出さず。

*** seismic-noise-tomography *************
https://github.com/bgoutorbe/seismic-noise-tomography
Rayleigh wave group velocities at periods of 6–23 s across Brazil from ambient noise tomography
・Python2.7でした。ライブラリがそろっていた 3.8 環境で動かし始めましたが、変換が面倒になり断念。

*** amb_noise_processing *****************
http://hestia.lgs.jussieu.fr/~boschil/downloads.html
Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations
・これもPython2.7でした。変換が print 文のみでしたので、3.8 にて実施。

この他、MATLAB を利用した GUI 対応コードもいくつかありました。さすがに購入しようとは思いませんでしたが、仕事で必要になれば選択肢に入るでしょう。

私の環境では、最後の amb_noise_processing のみ分散曲線まで確認できました。詳細を追いましょう。

************************************
20210516追加
「Two-receiver measurements of phase velocity: cross-validationof ambient-noise and earthquake-based observations」より、passive 手法で共通と思われる部分。

• Microseisms are most energetic in the period band between 5 and 30 s. At the shorter period end of this range, surface waves from teleseismic earthquakes are difficult to observe. Local to regional seismicity is needed to obtain phase velocities from earthquake data in the frequency range between about 3 and 10 s (Endrun et al. 2004).
• 1年のデータが必要（１年以上は不要）
• Records are filtered in the frequency band of interest (2–200s). サンプリング周波数1Hz以下⇒容量削減、計算負荷軽減。
• The data is then cut into overlapping time windows (e.g. half an hour with 50 per cent overlap), cosine-tapered at both ends.

 

C言語の基礎

GW 中に、C言語を眺めておこうと考えていました。

英語同様に、辞書があれば読めるものの全く書けない、という状況です。必要に迫られることが少なく、積極的に利用する機会はありませんでした。が、書ける能力があれば楽だったかも、ということは数回ありました。将来？に備えた弱点補強です。

取り急ぎ作りたいモノもないので、いくつかのサイトを見て１週間程度で完了できそうな以下を実施。

一週間で身につくC言語の基本
http://c-lang.sevendays-study.com/index.html

基礎編はほぼ復習になりましたが、応用編は知らないことばかり。あわせてメモしておきましょう。

基礎編
・void とは、関数の戻り値が無いことを意味する。
・Visaul Studio では、セキュリティの関係で scanf を使えない（scanf_s　ならOK）。
 プロジェクトを右クリック → プロパティ → c/c++ → SDL チェックを「いいえ」に変更
　※画面に直接入力する仕様は少ない（が、たまにあると忘れる）のでメモ。
・７日目、ファイル分割の御作法。
　#ifndef、#define、#endif マクロの意味、２重インクルードの防止、global 変数に対する extern修飾子。

応用編
２日目　変数のアドレス
・通常の変数：int a    値a    アドレス&a
・ポインタ変数：int* p    値*p    アドレスp
    p = &a:    pに、aのアドレスを代入    *pはaと同じものになる
    *p = 300:    *pに300を代入
3日目　ポインタと配列
・ポインタへのアドレス代入：p1 = &ar1[0]
    配列変数    アドレス    ポインタ    ポインタ変数の値
    ar1[0]    &ar1[0]    p1        *p1
    ar1[1]    &ar1[1]    p1+1    *(p1+1)
・ポインタ変数ｐ1に配列d[3] = { 0.2 , 0.4 , 0.6 }を代入：配列変数dの先頭アドレス&d[0]の値を代入するのと同じ処理
・動的メモリ確保 malloc()～解放 free()を使う場合、stdlib.hをインクルード
5日目　構造体
・構造体の成分の変数「メンバ」：(構造体変数名).(メンバ)
・データ渡し⇒スタック領域の圧迫、データコピー（無駄な処理）発生。⇒ポインタ渡しが通常

初心者がまず躓くといわれているポインタ。shallow・deep copy の関係と似ているようです。概念は理解できましたが、慣れるまでに時間がかかるでしょう。

丁寧な御説明で分かりやすく、2日で終わりました。感謝。
今回はココまで。

粒子法入門 その2

昨日読み終わってしまい、時間があったのでソースファイルを動かしてみました。

土砂移動のTESTとしてソースを触るにしても、オリジナルを動かせないと始まりません。
C++初心者ですが、ま、この程度は簡単でしょう、と簡単に考えていました。が、半日かかりました。

主な原因は、コンパイラーの違い。Linux ＋ gcc で動くように書かれていたのですが、私の環境はWin8.1 + VS_Express の VC++2013。そのままでは走らないことに、手を動かしてはじめて気づきました。といっても、押さえどころは2箇所です。


1　fopen などは 2013 でそのままでは使えません。今回はセキュリティー上お奨めの fopen_s などに書き換えました。

//fp = fopen(filename, "r");
//fscanf(fp,"%d",&np);

errno_t err;
err = fopen_s(&fp, filename, "r");
if (err == 0)
{
fscanf_s(fp, "%d", &np);

このあたりは多くのブログに載っていますね。戸惑いやすい点だったのでしょう。


2　時間計測。そのままでは動かなかったので、ヘルプを見ながら書き換えてみました。間違っているかもしれませんが、時間計測は個人的に使わないので（いつもの計算は数日かかるので出力ファイルのタイムスタンプで見ています）、動けば良い程度に。

/*#include <sys/time.h>
double get_dtime(void){
struct timeval tv;
gettimeofday(&tv, NULL);
return ((double)(tv.tv_sec) + (double)(tv.tv_usec) * 0.000001);
}*/
#include <sys/timeb.h>
#include <sys/types.h>
double get_dtime(){
struct _timeb tv;
_ftime64_s(&tv);
return ((double)(tv.time) + (double)(tv.millitm) / 1000.0);
}


ここを超えると、あとは簡単。

粒子配置版を実行し、できたファイルをOpenMP 版で計算し、、ParaView 用に変換、可視化する流れです。以前と同じ、ダムブレイクです。SPHysics の時は粒子数いくつだったでしょうか？今回は同じ環境ですが、粒子数1.9万、計算時間は数秒でした。
http://phreeqc.blogspot.jp/2014/01/paraview.html

とりあえずは、動きました。入門できたでしょうか？

OpenCV + VC++ の設定

OpenCV 2.2を利用するため、Visual C++ 2010 Express を インストールしました。

基本的には、VC++も、VB や VF と似たような構成です。OpenCV を利用するだけなら、それほど深い知識は必要なさそうです。

空のWin32コンソールアプリを作成し、コードを書き込みます。このままでは OpenCV へのリンクが設定されていませんので、以下のパスを追加します。

「プロジェクト」-「プロパティ」-「構成」-「全ての構成」
「プロジェクト」-「プロパティ」-「C/C++」-「全般」-「追加のインクルードディレクトリ」
    C:\Program Files\OpenCV\include　＜＜＜これが必要！
    C:\Program Files\OpenCV\include\opencv
    C:\Program Files\OpenCV\include\opencv2
その他必要に応じて。

「プロジェクト」-「プロパティ」-「リンカー」-「全般」-「追加のライブラリディレクトリ」
    C:\Program Files\OpenCV\lib

ライブラリはソース内で指定しました。
    #pragma comment(lib,"C:\\program files\\OpenCV\\lib\\opencv_core220d.lib")
    #pragma comment(lib,"C:\\program files\\OpenCV\\lib\\opencv_imgproc220d.lib")
    #pragma comment(lib,"C:\\program files\\OpenCV\\lib\\opencv_highgui220d.lib")
    #pragma comment(lib,"C:\\program files\\OpenCV\\lib\\opencv_objdetect220d.lib")
デバッグの場合はdのついたもの、リリースではdのないファイルが必要です（#ifdef 文で選択しま
す）。その他必要に応じて。

これでOpenCVが使えるようになりました。
optical flow のサンプルコードもビルド可能となり、２つの画像を比較して差分をベクトル表示できるようになりました。

------------------------------------------------------
2/1追記
ツール-設定-上級者用の設定
表示-プロパティ―マネージャーより debug - Microsoft.Cpp.＿＿＿.user をダブルクリックで、VC++ ディレクトリへの入力ができました。