物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

Python scipyのcurve_fit で導出したパラメータと標準誤差をCERN ROOTの結果と比較

まずはCERN ROOT + C++ で実装する。

お手本通り、平均値0、標準偏差1、ガウス分布(正規分布)に沿う乱数を10000個作り、ROOTのヒストグラムに詰め、 TF1 の ガウス分布 gaus でフィッティングした。オプション等は何もつけていない。\chi^{2}も普通に結果を引用した。

正規分布の関数は、


f(x)=const. \exp \! \left( -\frac{(x-mean)^2}{sigma^2} \right)

である。

Constant 402.531 +/- 4.95816
Mean     -0.00397889 +/- 0.00989405
Sigma    0.983115 +/- 0.00704207
Chi2     83.6418

次に、全く同じデータを使うために出力しておいたテキストデータを使って、 scipy.optimize.curve_fit でフィッティングさせる。 テキストデータはC++で発生させることが可能だが、C++のコンパイル環境がない人のためにファイルも提供する。

ガウス分布のフィッティング用のサンプルデータ Sample data

ガウス分布の式は自分で作成し与えた。また、\chi^{2}scipy.stats.chisquare を使って計算した。ここでポイントは、誤差を適切に与えることである。ビンあたりのエントリー数Nに対して、誤差\sqrt{N}の配列を curve_fitsigma に与えるとよい。

Chi2 83.641822
            Estimate  Std. error
Constant  402.530641    4.963109
Mean       -0.003979    0.009894
Sigma       0.983115    0.007064

推定値は完全に一致。誤差は若干異なるがなぜだろう。愛嬌でしょうか。誤差の3桁目までは一致しているので、実用上は問題ないだろう。 scipy.optimize.curve_fit においても、誤差を適切に与えることでROOTの結果と推定値は完全に一致することが分かった。

p値(p-value)を求めるときは、カイ二乗をc、自由度数をndfとすると、ROOTだとTMath::Prob(c, ndf)、Pythonだと

from scipy.stats import chi2
chi2.sf(c, ndf)

とする。今回の例だと TMath::Prob(83.641822, 69)=chi2.sf(83.641822, 69)=0.11055024 である。

以下ソースコード。

C++ と ROOT版 C++ と ROOT5版 ガウス分布をフィッティング ($1781181) · Snippets · GitLab

Python3版 Python3版 ガウス分布をフィッティング ($1781185) · Snippets · GitLab

ウェブの情報をいっぱい参考にしたが、その一部。

https://omedstu.jimdo.com/2018/07/16/python%E3%81%A7gaussian-fitting/

scipy.stats.chisquare — SciPy v1.11.1 Manual

講演会 川村静児氏「重力波:アインシュタインの奏でる宇宙からのメロディー」質疑応答 @2018年 ぎふサイエンスフェスティバル

会場内でメモしたため間違っているところがあるかもしれません。

Q. インフレーション理論で物理現象(特殊相対性理論)は成立しているのか?
A. 場が広がってるだけなので成立できる。
Q. インフレーションの音はシミュレーションはされているのか?
A. されているが、まだ言わない。
Q. 重さがなくなったのは重力子の質量になったと考えられないのか?
A. 重力波の媒介粒子は質量はないと考えられているので、ほぼ全てがエネルギーになったと考えるのが妥当。
Q. ブラックホール連星では他の情報は出てこないのか?
A. ブラックホール連星の合体では出てこないとされている。中性子星連星の合体では他の情報も出てくる。
Q. 光と重力波は何が違うのか?
A. 媒介する力が違う。また、電磁波はベクトル波、重力波はテンソル波。
Q. 重力子と重力波は同じなのか?
A. 重力波も粒であり波であるから、重力子と考えても良いと思う。

注: インフレーションで発生する重力波はsub-Hz帯なので、実は人間の耳には聞こえないんです。中性子連星やブラックホール連星での合体で発生する重力波は たまたま 人の可聴領域内だっただけなのです。

http://www.city.gifu.lg.jp/17404.htm

川村さん、準備と講演、お疲れさまでした。

updatestar.com について

このサイトから、2つほどファイルをダウンロードしたところ、一つはMicrosoft Defender君がきっちりガードしてくれた。

別のファイルをダウンロードしたところ、Microsoft Defender は動かなかったが、SHA-256を計算しVirusTotalで検索してみたところ見事引っかかった。

VirusTotal (←これは真っ当なサイトよ!)

たぶん、いえ、間違いなく、updatestarはダメな奴です。

検索用: コンピューターウイルス・アンチウイルスソフト

OpenCV 2系でGPUで膨張処理をさせる方法

OpenCV 2系でGPUで膨張処理をさせる方法で詰まったので書いておく。最後に検証に使った全コードがあるのでどうぞ。

CPUで膨張処理をさせるとき、さくっと書けば次のようになる。

cv::Mat src, dst;
int width = 256;
int height = 256;
src = cv::Mat::eye(cv::Size(width, height), CV_8U);
cv::dilate(src, dst, cv::Mat::ones(3, 3, CV_8U));
std::cout << cv::countNonZero(dst) << " cpu dilate w/o rect" << std::endl;

この結果は、 256+255*2+254*2=1274 となる。 たとえば、周辺1ピクセルを膨張処理させなかった場合は、

cv::Rect rect = cv::Rect(1, 1, width - 2, height - 2);
cv::dilate(src(rect), dst(rect), cv::Mat::ones(3, 3, CV_8U));

この結果は、 256+255*2+254*2-5*2=1264 となる。

GPU版では、領域外チェックが甘いため、全ての領域を処理させたとき、領域外のメモリを含めて計算してしまうバグ(仕様)がある。 それを検証するために、 cv::gpu::dilate と、 cv::gpu::getMaxFilter_GPU と、cv::gpu::createBoxFilter_GPU を比較してみた。 それぞれ、全領域で処理させた場合 w/o rect と、 rect で周辺1ピクセルを無視した場合 with rect を行った。 結果は次のとおり。

処理方法 countNonZeroの結果
gpu dilate w/o rect 1264
gpu dilate with rect 1249
gpu max w/o rect 1278
gpu max with rect 1264
gpu box w/o rect 1259
gpu box with rect 1264

w/o rect の期待値は1274なので、いずれも正しい答えを返していない。一方で、 with rectcv::gpu::getMaxFilter_GPU と、cv::gpu::createBoxFilter_GPU で正しい答えを返している。 この2つ以外は、おそらく環境依存で値が変わりうると考えられる。ここで、BoxFilterを使うときは、下記の検証用コードでもそそしているとおり、要素を Expansion の二乗より大きい値にしておく必要がある。さらに、最後に適切な値に閾値処理する必要がある。これらの手間を考えると、 cv::gpu::getMaxFilter_GPUwith rect で使うのが最もシンプルな方法であろう。

なかなか難解なライブラリである。

さて、検証用全コード

#include <opencv2/opencv.hpp>
#include <opencv2/gpu/gpu.hpp>

using namespace cv;

int main() {

    cv::gpu::setDevice(0);

    int gpudeviceid = cv::gpu::getDevice();
    cv::gpu::DeviceInfo dev(gpudeviceid);
    string gpu_name = dev.name();

    cv::gpu::Stream stream;
    int Expansion = 3;

    Ptr<cv::gpu::BaseFilter_GPU> fb;
    fb = cv::gpu::getMaxFilter_GPU(CV_8UC1, CV_8UC1, Size(Expansion, Expansion));

    cv::Ptr<gpu::FilterEngine_GPU> filter;
    filter = gpu::createBoxFilter_GPU(CV_8U, CV_8U, cv::Size(Expansion, Expansion));

    cv::gpu::GpuMat gsrc, gdst;
    cv::Mat src, dst;

    int width = 256;
    int height = 256;
    src = cv::Mat::eye(cv::Size(width, height), CV_8U) * 9;
    gsrc.upload(src);

    cv::Rect rect = cv::Rect(1, 1, width - 2, height - 2);

    {
        cv::dilate(src, dst, cv::Mat::ones(3, 3, CV_8U));
        std::cout << cv::countNonZero(dst) << " cpu dilate w/o rect" << std::endl;
    }

    {
        dst = cv::Scalar(0);
        cv::dilate(src(rect), dst(rect), cv::Mat::ones(3, 3, CV_8U));
        std::cout << cv::countNonZero(dst) << " cpu dilate with rect" << std::endl;
    }

    src = cv::Scalar(0);
    {
        gdst.upload(src);
        cv::gpu::dilate(gsrc, gdst, cv::Mat::ones(3, 3, CV_8U));
        std::cout << cv::countNonZero(dst) << " gpu dilate w/o rect" << std::endl;
    }

    {
        gdst.upload(src);
        cv::gpu::dilate(gsrc(rect), gdst(rect), cv::Mat::ones(3, 3, CV_8U));
        std::cout << cv::gpu::countNonZero(gdst) << " gpu dilate with rect" << std::endl;
    }

    {
        gdst.upload(src);
        (*fb)(gsrc, gdst);
        std::cout << cv::gpu::countNonZero(gdst) << " gpu max w/o rect" << std::endl;
    }

    {
        gdst.upload(src);
        (*fb)(gsrc(rect), gdst(rect)); //おすすめの方法
        std::cout << cv::gpu::countNonZero(gdst) << " gpu max with rect" << std::endl;
    }

    {
        gdst.upload(src);
        filter->apply(gsrc, gdst);
        std::cout << cv::gpu::countNonZero(gdst) << " gpu box w/o rect" << std::endl;
    }

    {
        gdst.upload(src);
        filter->apply(gsrc, gdst, rect);
        std::cout << cv::gpu::countNonZero(gdst) << " gpu box with rect" << std::endl;
    }

    fb.release(); // no error

    stream.waitForCompletion();
    std::cout << "Free memory rate = " << dev.freeMemory()*100.0 / dev.totalMemory() << " [%] at Device = " << gpudeviceid << " at first" << std::endl;
    return 0;
}

WSL (Ubuntu系)における便利なコマンド一覧

Windows Subsystem on Linux

自分用のメモとして随時更新

インストールパッケージをアップデート。これをしないと何もできない。

sudo apt update
sudo apt upgrade

Anacondaをインストールする

sudo sh Anaconda3-2018.12-Linux-x86_64.sh

ビルドツールを一括でインストール。gccやmakeなどがインストールされる。

sudo apt install build-essential

ファイルをダウンロード

wget https://ファイル名

debファイルのインストール

sudo dpkg -i code_1.32.3-1552606978_amd64.deb

なんか壊れた時

sudo apt --fix-broken install

gccのバージョン

gcc -v

bashを実行する時に読み込まれるもの。bash_profileとbashrcは何がどう違うの?

pico ~/.bash_profile
source ~/.bash_profile
pico ~/.bashrc
source ~/.bashrc

ROOTのインストール

export ROOTSYS=/(インストールしたフォルダ)/root_v6.16.00.Linux-ubuntu18-x86_64-gcc7.3
export PATH=$PATH:$ROOTSYS/bin
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$ROOTSYS/lib
# echo $ROOTSYS

" スプラッシュを表示させないには以下を記述

alias root="root -l"

.rootlogon.C の例

https://gitlab.com/yoshimoto/root_macros/commits/master/.rootlogon.C

~/.rootrc の例。GuiのFontサイズを変更できる。

Gui.DefaultFont:            -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-1
Gui.MenuFont:               -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-1
Gui.MenuHiFont:             -*-helvetica-bold-r-*-*-12-*-*-*-*-*-iso8859-1
Gui.DocFixedFont:           -*-courier-medium-r-*-*-12-*-*-*-*-*-iso8859-1
Gui.DocPropFont:            -*-helvetica-medium-r-*-*-12-*-*-*-*-*-iso8859-1
Gui.IconFont:               -*-helvetica-medium-r-*-*-10-*-*-*-*-*-iso8859-1
Gui.StatusFont:             -*-helvetica-medium-r-*-*-10-*-*-*-*-*-iso8859-1

ROOTをグラフィックスなしで起動する

ROOT -b

GUI用の

sudo apt install xfce4-terminal
sudo apt install xfce4

# .bash_profile または .bashrcに
export DISPLAY=:0.0
# high DPIのときは
export GDK_SCALE=2
export GDK_DPI_SCALE=2

# GUIのテストには
xeyes

VcXsrv.exe(おそらく C:\Program Files\VcXsrv\ あたりにある) のプロパティから 高DPI設定の変更 → 高いDPIスケールの動作を上書きします。 のチェックを入れる

[備忘録]高DPI環境での WSL+VcXsrv の設定 - Qiita

ただし、一部アプリでは無効なので、あまりご利益はないかもしれない。(例: ROOT)

WSLにSSH接続する

https://qiita.com/ezmscrap/items/30eaf9531e240c992cf1

sudo apt install openssh-server
sudo nano /etc/ssh/sshd_config
# PasswordAuthentication no ->
# PasswordAuthentication yes
sudo service ssh restart

gdbを起動できませんでした。システムにgdbが見つからないため、インストールが必要です。システムのパッケージマネージャーを使用してインストールしてください。

と出るので、インストールする

sudo apt install gdb

Poco C++ Librariesをインストールする

https://pocoproject.org/docs/00200-GettingStarted.html#9

Javaのインストール

sudo add-apt-repository ppa:webupd8team/java

sudo apt update
sudo apt install oracle-java8-installer

NVIDIAのドライバをインストールする

https://packages.ubuntu.com/ja/bionic/nvidia-384

384を入れる時は

sudo apt update
sudo apt install nvidia-384
nvidia-smi

最新のを入れる時は

sudo add-apt-repository ppa:graphics-drivers/ppa
sudo apt update
sudo apt install nvidia-390
sudo apt install nvidia-410

両立できないので、どちらかだけよ。

ホームディレクトリにユーザーがあるか確認

ls -l /home

ユーザーのパスワードがあるかで確認

cat /etc/passwd

ユーザーの追加

sudo adduser newusername

ユーザーの削除

sudo userdel -r deletableusername

ウェブページの一括ダウンロード

wget -r -l 10 https://URL/
#-r --recursive 再帰ダウンロードを行う(URLで指定したファイルからのリンクもダウンロードする。デフォルトは5階層まで:実行例3を参照)
#-l 数 --level=数 指定した回数分リンクをたどる(「-r」オプションは「-l 5」相当)