読者です 読者をやめる 読者になる 読者になる

Physics-Station 2

旧 http://physics-station.blogspot.jp/ から当はてなブログに移行しました。間違ってるところがあればコメントください。記述の正確性は保証しません。

OpenCVのgpu::countNonZeroをgpu::Streamで高速化する

cv::gpu::countNonZerogpu::Streamで高速化したいという需要があったので、その解決策を書いておく。

opencvの2系のcountNonZeroはstreamを使うことが出来ない。そのため、非同期処理を行い高速化する際のボトルネックになる。cv::gpu::reduceは2次元画像を1次元に射影するための関数でstreamが使えるので、CV_REDUCE_SUMを使って射影した後、軽量化された1次元画像を非同期ダウンロードし、CPUで積算処理を行えばよい。行または列の合計値が255を超えた場合は完全に一致しないものの、ラフな値で良い場合はこの方法が適切である。

以下はサンプルコードでは、コアの部分だけ記述した。streamやgvec等を使いまわさないと初期化に時間がかかるので注意が必要。

cv::gpu::GpuMat gmat = cv::gpu::GpuMat(cv::Size(1920, 1080), CV_8UC1);
cout << cv::gpu::countNonZero(gmat) << endl; //streamが使えない

cv::gpu::GpuMat gvec;
cv::gpu::Stream stream;
cv::gpu::reduce(gmat, gvec, 1, CV_REDUCE_SUM, -1, stream); //行ごとに積算する場合
cv::Mat vec;
stream.enqueueDownload(gvec, vec);
cout << cv::sum(vec)[0] << endl;

C# でWindows 上のプログラムのリソースを監視する

C#で以下のリソースの監視をしたいという需要があったのでサンプルコードを公開する。

  • 使用メモリ量
  • ハンドル数
  • GDI オブジェクト数
  • User オブジェクト数
using System;
using System.Diagnostics;
using System.Runtime.InteropServices;

namespace ConsoleApplication1
{
    class Program
    {
        [DllImport("user32.dll")]
        private static extern uint GetGuiResources(IntPtr hProcess, uint uiFlags);
        const uint GR_GDIOBJECTS = 0;
        const uint GR_USEROBJECTS = 1;
        /// <summary>
        /// GDI オブジェクトの数を返す
        /// </summary>
        /// <returns></returns>
        public static uint GetGDIObjects()
        {
            return GetGuiResources(Process.GetCurrentProcess().Handle, GR_GDIOBJECTS);
        }
        /// <summary>
        /// User オブジェクトの数を返す
        /// </summary>
        /// <returns></returns>
        public static uint GetUserObjects()
        {
            return GetGuiResources(Process.GetCurrentProcess().Handle, GR_USEROBJECTS);
        }
        static void DisplayMemory()
        {
            Console.WriteLine("Total memory  : {0:###,###,###,##0} bytes", GC.GetTotalMemory(false));
            Console.WriteLine("Private bytes   {0:###,###,###,##0} bytes", Process.GetCurrentProcess().PrivateMemorySize64); //プライベート メモリの量
            Console.WriteLine("Handle   count: {0}", Process.GetCurrentProcess().HandleCount); // ハンドル数
            Console.WriteLine("GDI  obj count: {0}", GetGDIObjects());
            Console.WriteLine("User obj count: {0}", GetUserObjects());
            Console.WriteLine();
        }
        static void Main(string[] args)
        {
            for (int i = 0; i < 100; i++)
            {
                using (var p = new Process())
                {
                    p.StartInfo.FileName = "ipconfig";
                    p.Start();
                    p.WaitForExit();
                }
            }
            DisplayMemory();
            GC.Collect(); //ガーベッジ コレクションを行う
            GC.WaitForPendingFinalizers(); //ガーベッジ コレクションが終わるまで待つ
            DisplayMemory();
        }
    }
}

標準出力とタスクマネージャーはこんな感じ

f:id:onsanai:20161004160344p:plain:w600

参照

.net - Finding Memory leaks in C# - Stack Overflow

GDI object count - C# / C Sharp

Windows 10のアップデート後の不本意な再起動を止める方法

f:id:onsanai:20161001184332p:plain:w600

Windows Updateを生かしたまま、アップデート後の不本意な再起動を止める方法は存在しない。Pro版でも、Enterprise版でも。

なので、Windows Updateのサービスを止めるしかない。という結論にたどり着いた。

諦めよう、Microsoftが改心するまで。

blog.bagooon.com

3Dの位置と角度を持った情報の2次元プロジェクションマップ

CERNが開発しているROOT (バージョン5)を用いた検出器内の飛跡を2次元プロジェクションに投影して描画する方法。 具体的には三次元情報を持つ飛跡を二次元に投影したい時に使う。

#include <vector>
#include <random>
#include <limits>
#include <TArrow.h>
#include <TGraph.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TPaveText.h>
class Track
{
public:
    double px, py, ax, ay;
};
int main()
{
    std::random_device rand;
    std::vector<Track> vtrack;

    double width = 10;
    double height = 10;

    TCanvas *c1 = new TCanvas("c1", "", 1000, 1000);
    
    for (int i = 0; i < 100; i++) {
        Track t;
        t.px = double(rand()) / std::numeric_limits<unsigned int>::max() * width;
        t.py = double(rand()) / std::numeric_limits<unsigned int>::max() * height;
        t.ax = double(rand()) / std::numeric_limits<unsigned int>::max() - 0.5;
        t.ay = double(rand()) / std::numeric_limits<unsigned int>::max() - 0.5;
        vtrack.emplace_back(t);
    }
    double extz = 2; //ファクター
    double arrowsize = 0.01;

    TGraph* dummy = new TGraph();
    dummy->SetPoint(0, 0, 0); //左下
    dummy->SetPoint(1, width , height ); //右下
    dummy->SetMarkerColor(0); //マーカー色は白に

    dummy->GetXaxis()->SetTitle("X");
    dummy->GetXaxis()->SetLabelSize(0.04);
    dummy->GetXaxis()->SetTitleSize(0.04);

    dummy->GetYaxis()->SetTitle("Y");
    dummy->GetYaxis()->SetLabelSize(0.04);
    dummy->GetYaxis()->SetTitleSize(0.04);

    dummy->SetTitle("Arrow");
    dummy->Draw("AP"); //軸とポイントを描画

    for (auto p : vtrack) {

        TArrow *arrow = new TArrow(p.px, p.py, p.px + (p.ax)*extz, p.py + (p.ay)*extz, arrowsize, ">");
        arrow->SetFillColor(kBlue);
        arrow->SetLineColor(kBlue);
        arrow->SetFillStyle(1001);
        arrow->SetLineWidth(1);
        arrow->SetLineStyle(1);
        arrow->Draw();
    }

    //リファレンスの矢印
    double ref_px = -width / 8, ref_py = -height / 8;
    {
        TArrow *arrow = new TArrow(ref_px, ref_py, ref_px + extz, ref_py, arrowsize, ">");
        arrow->SetFillColor(kRed);
        arrow->SetLineColor(kRed);
        arrow->SetFillStyle(1001);
        arrow->SetLineWidth(1);
        arrow->SetLineStyle(1);
        arrow->Draw();
    }
    {
        TArrow *arrow = new TArrow(ref_px, ref_py, ref_px, ref_py + extz, arrowsize, ">");
        arrow->SetFillColor(kRed);
        arrow->SetLineColor(kRed);
        arrow->SetFillStyle(1001);
        arrow->SetLineWidth(1);
        arrow->SetLineStyle(1);
        arrow->Draw();
    }
    //リファレンスの文字
    double text_px = -width / 10, text_py = -height / 10;
    {
        TPaveText *pt = new TPaveText(text_px, text_py, 0, text_py + height / 30, "br");
        pt->SetFillColor(0);
        pt->SetLineColor(0);
        pt->SetShadowColor(0);
        TText *text = pt->AddText("1.0 [rad]");
        pt->Draw();
    }

    c1->SaveAs("Arrow.pdf");

    return 0;
}

作った画像

f:id:onsanai:20160919174322p:plain

2個の数列(std::vector<int>とか)から重複とかを探す

2個の数列から重複等を探す方法。

続きは 2個の自作クラスの配列(std::vector<MyClass>とか)から重複とかを探す - Physics-Station2

#include <algorithm>に便利な関数が用意されている。 ソートして、重複を削除した後、set_intersection, set_union, set_differenceを使ってみよう。

参照

Tech Tips: 知ってると便利なSTL(10) set_intersection, set_union, set_difference

乱数発生器についてのコードを修正しました(2016/11/30)

list1 = {5,7,4,2,9}
list2 = {4,3,8,3,2}
↓ソート
list1 = {2,4,5,7,9}
list2 = {2,3,3,4,8}
↓重複を削除
list1 = {2,4,5,7,9}
list2 = {2,3,4,8}
↓比較
積集合 = {2,4}
和集合 = {2,3,4,5,7,8,9}
差集合1 = {5,7,9}
差集合2 = {3,8}
#include <random>
#include <vector>
#include <algorithm>
#include <iostream>
#include <iterator>

void main()
{
    const int Ntrk = 10000000;

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> dis(0, Ntrk - 1);
    std::vector<unsigned int> list1, list2;

    for (int i = 0; i < Ntrk; ++i) {
        list1.emplace_back(dis(gen)); //乱数を作成
        list2.emplace_back(dis(gen)); //乱数を作成
    }

    std::sort(list1.begin(), list1.end()); //ソート
    std::sort(list2.begin(), list2.end()); //ソート

    std::cout << list1.size() << " -> ";
    list1.erase(std::unique(list1.begin(), list1.end()), list1.end()); //重複を削除
    std::cout << list1.size() << std::endl;

    std::cout << list2.size() << " -> ";
    list2.erase(std::unique(list2.begin(), list2.end()), list2.end()); //重複を削除
    std::cout << list2.size() << std::endl;

    std::vector<unsigned int> list_inter; //list1とlist2の両方にある値 (積集合)
    std::set_intersection(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_inter));

    std::vector<unsigned int> list_union; //list1かlist2のどちらかに存在する値(和集合)
    std::set_union(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_union));

    std::vector<unsigned int> list_dif1; //list1にあってlist2にない値(差集合)
    std::set_difference(list1.begin(), list1.end(), list2.begin(), list2.end(), std::back_inserter(list_dif1));

    std::vector<unsigned int> list_dif2; //list2にあってlist1にない値(差集合)
    std::set_difference(list2.begin(), list2.end(), list1.begin(), list1.end(), std::back_inserter(list_dif2));

    std::cout << "dif1 + dif2 + inter = union" << std::endl;
    printf_s("%u + %u + %u = %u\n", list_dif1.size(), list_dif2.size(), list_inter.size(), list_union.size());
    ::system("pause");
}

δtanθ - tanθグラフをδθ - tanθグラフに変換する

より正しい記述へ

phst.hateblo.jp

tanθを0.1刻みで、δtanθの値が一律で0.1のときにδθを計算するコードの例。

#include <vector>

void main()
{
    std::vector<double> x, y; //tanθ δtanθ
    for (int i = 0; i < 20; i++) {
        x.emplace_back(i * 0.1 + 0.05);
        y.emplace_back(0.1);
    }

    for (int i = 0; i < 20; i++) {
        double dtheta = atan(x.at(i) + y.at(i) / 2) - atan(x.at(i) - y.at(i) / 2); //δθを計算
        printf("%8.5f %8.5f %8.5f\n", x.at(i), y.at(i), dtheta);
    }
}

f:id:onsanai:20160918044320p:plain

OpenMPでfor文を高速化する

openmpでfor文を高速化してみよう。OpenMPを有効にするには、Visual Studioのプロジェクトのプロパティページを開いて、C/C++の言語のOpenMPのサポートをはい(/openmp)にする必要がある。

例として1から50000までの数が素数(prime)かどうかを調べてみよう。素数の調べ方は2以上候補値未満の数で割った余りを調べる愚直な方法で実装した。

#include <omp.h>
#include <vector>
#include <iostream>
#include <chrono>
#include <algorithm>

int main()
{
    auto *mylock = new omp_lock_t;
    omp_init_lock(mylock); //lockを初期化

    std::vector<int> result; //素数の結果

    auto start = std::chrono::system_clock::now(); //開始時間

#pragma omp parallel for num_threads(12)
    for (int i = 0; i < 100; i++) {
        std::vector<int> list; //素数の中間ファイル
        
        for (int j = 0; j < 5000; j++) {
            int candidate = i + j * 100 + 1;
            if (candidate < 2) { continue; }

            bool prime_flag = true;
            for (int k = 2; k < candidate; k++) {
                if ((candidate%k) == 0) { //2以上候補値未満の数で割った余り
                    prime_flag = false; //素数ではない
                    break;
                }
            }
            if (prime_flag) {
                list.emplace_back(candidate);
            }
        }
        omp_set_lock(mylock); //std::vector<T>の変更操作をするためlock
        for (auto p = list.begin(); p != list.end(); p++) {
            result.emplace_back(*p);
        }
        omp_unset_lock(mylock); //lockを解除
    }
    auto end = std::chrono::system_clock::now(); //終了時間

    omp_destroy_lock(mylock); //lock破棄

    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " [msec]" << std::endl;

    std::cout << "max prime is " << *std::max_element(result.begin(), result.end()) << std::endl; //調査した中で最大の素数を出力
    ::system("pause");
}

当環境で 37秒かかったものが12スレッド並列で行うと6.5秒に短縮された。オーバーヘッドがあるため、単純に1/12にはならない。 std::vector<T>の変更操作はスレッドセーフではないため、omp_lock_tを使ってlockとlock解除を行っている。

素数の調べ方は最適化しておらず、世の中にはもっと高速な方法が沢山ある。一般的に、openmpで高速化するのはアルゴリズムを改善した後のほうが良いだろう。

スレッド数を指定しない場合は以下のように記述することも可能。この場合は全スレッドを使って計算が行われる。

#pragma omp parallel for

Visual Studio 2013でstd::min std::maxが使えなくなった時の対処

std::max - cppreference.com

Visual Studio 2013以降、Visual Studio 2015も完全にalgorithmに移動したらしいので以下の記述が必要。

#include <algorithm> //必要
int main()
{
    int i = std::max(1,2);
    int j = std::min(3,4);
}

ちなみにVisual Studio 2012以下は以下のコードで動く。

#include <xutility> //非推奨
int main()
{
    int i = std::max(1, 2);
    int j = std::min(3, 4);
}

VC10, VC11ではalgorithm 経由で xutilityに記述、VC12, VC14ではalgorithm に直接記述されるようになったようだ。VC12, VC14では、別のincludeファイルからxutilityが参照されている場合は、明示的にalgorithmを追加してやる必要がある。

メモリリークしているソースコードの行を特定する

#include <cstdlib>
#include <new>
#include <memory>

#include <crtdbg.h>

#define _CRTDBG_MAP_ALLOC

#define new ::new(_NORMAL_BLOCK, __FILE__, __LINE__)

void main()
{
    int *i = new int;
    _CrtDumpMemoryLeaks();
    return;
}

これをデバックモードで実行すると

Detected memory leaks!
Dumping objects ->
c:\users\documents\source.cpp(13) : {170} normal block at 0x00000200AD5747B0, 4 bytes long.
 Data: <    > CD CD CD CD 
Object dump complete.

ちゃんとint *i = new int; である13行目でエラーが出ていることが分かる。

参照

プログラマーの友 第八報:メモリリークと crtdbg.h

catchしなかった例外が発生した際の ***は動作を停止しました「WerFault」を表示させない

f:id:onsanai:20160912070612p:plain

問題が発生したため、プログラムが正しく動作しなくなりました。プログラムは閉じられ、解決策がある場合はWindowsから通知されます。

を表示させないプログラムは以下のとおり。

#include <exception>
#include <Windows.h>
void main()
{
    SetErrorMode(SEM_NOGPFAULTERRORBOX);
    throw std::exception(); //例外を投げる
}

参照

SetErrorMode function (Windows)

逆引きWIN32API: 一般保護例外ダイアログを出さなくする方法 - seclan のほえほえルーム

[VC]子プロセスが落ちたのを検出する。(6) うずまき の なんとなくでいいのかも?

マルチスレッドにおける例外処理の受け渡し (VC++)

別スレッド中の例外を、本スレッドに渡す方法のメモ。動作確認はVisual Studio 2013 C++で行っている。

詳しい説明は スレッド間の例外転送 を参照してね。

#include <thread>
#include <iostream>

//マルチスレッド用の関数
void f1(std::exception_ptr &eptr) {
    try {
        throw std::exception("test");
    }
    catch (...) { //すべての例外を受ける
        eptr = std::current_exception();
    }
}

int main() {
    std::exception_ptr eptr;
    std::thread t1(f1, std::ref(eptr)); //スレッドを作成し参照でeptrを渡す
    t1.join(); //終了待ち

    try {
        if (eptr) {
            std::rethrow_exception(eptr); //例外を投げなおす
        }
    }
    catch (std::exception &ex) {
        std::cout << ex.what() << std::endl;
    }
    ::system("pause"); //デバッグ用のpause
}

WiMAX 2+は下り速度4.5Mbpsしか出ない

タイトルで何もかも終わったんだが、下り速度を測定してみた。まずは一番重要な速度の一つ、アニメのHD画質が見れるかどうか。

f:id:onsanai:20160904144227p:plain

確かにアニメのHD画質は見れるようだ。4.5Mbpsしか使っていないので、そもそも大したことはないんだが。

Re:ゼロから始める異世界生活のアニメ見放題 | dアニメストア

 

真面目に下り速度を測定してみる。

f:id:onsanai:20160904144535p:plain

6Mbpsか。これ、測定の時だけブーストしてないか?って疑いがあるので、信用しない方がいいと思う。

 

とにかく、速度制限に入ると本当に下り速度は4.5Mbpsになるから、WiMAX 2+を家庭用固定回線代わりに使うのはやめよう。遅すぎる。

 

次回、応答速度を測定予定

高速・安全・格安100TB ストレージを構築する

高速・安全・格安100TBストレージを構築してみよう。

安全なデータ領域を確保するためには、最低でもRAID1が必須だろう。

RAID 10を構築したいが、100TBもの領域をRAID0で構築しようとすると実データ領域で4倍必要になる。100TBの領域のために400TBの領域を確保するのは馬鹿らしい。

RAID 5では、大量のデータを高速で書き込む時に急激に速度が遅くなる問題が発生したので、ここでは除外する。

 

RAIDカードの選定だが、市販のRAIDカードでそこそこのクオリティなのがこの製品。

HighPointのページに詳細は書いてある。6 Gb/sで24ポート。バスインターフェースはPCIe 2.0 x 16。ポートはSFF-8087 (Mini-SAS) x6なのでこのタイプのケーブルが6本必要である。

HighPoint Global website

 

HDDはどうしよう。NAS対応縛りは必須だよね。

 

WD Red。以前は6TBまでだったが、執筆時点で8TBまで出たので容積あたりの容量がだいぶ増えた。35,725円(執筆時)。安いのはいいが、最近のHDD故障率レポートを見ると、安かろう悪かろうになりつつある。WD頑張れ。

Seagate IronWolf NAS。使ったことはないが10TBまで出ている。54,800円(執筆時)。今まで安かろう悪かろうのイメージだった海門ことSeagateの故障率が減っているのが特徴。

 WD傘下のHGST UltrastarAmazon価格で71,188円(執筆時)。故障率は非常に低いらしいが、高い。

 

ここではSeagateの10TBとWDの8TBを使ってみよう。

 

Seagateの10TB 54,800円を20台買って、110万円、RAID 1を考慮せずに200TB。

WDの8TB 35,725円を24台買って、86万円、RAID 1を考慮せずに192TB。

 

格安でいくならWDだけれど、故障率で交換のコストを考えるとSeagateの方がいいのかなあ。サーバー本体とRAIDカードを含めても130万ほどで収まるので、本当に安くなったものだ。

 

リンクはアフィを含むので注意。

10 Gbpsのテスト

10 Gbpsを導入したのでテストしてみた。

使ったのはIntelイーサネットサーバーアダプター X520-2とNETGEARのスイッチングハブGS752TXS。

www.intel.com

GS752TXS | Product | Support | NETGEAR

 

NETGEARのCuのSFP+のケーブル

とてもいいサーバー用パソコン

RAMディスク経由でCrystalDiskMark

以上。

 

さて、結果。

 

f:id:onsanai:20160819230537p:plain

 

CrystalDiskMarkを2プロセス同時に走らせて、左半分が2プロセスともRead、右半分が2プロセスともWrite、真ん中は片方がReadでもう片方がWrite。

 

送受信ともほぼ10 Gbpsを達成。うちもバックボーンが10 Gbpsになりました。気軽に10Gbpsを導入できる時代が来て感動。

 

次は、10GBASE-T対応のカテゴリ7のLANケーブルを使ってテストしてみる予定。

ELSA GeForce GTX 780 S.A.C を分解

 

 ELSA GeForce GTX780 S.A.Cを分解してみた

www.elsa-jp.co.jp

 

とても簡単に分解できた。全てのパーツを外すとこの通り。本体基盤の下がヒートスプレッダー、ヒートシンク、ファン・・・説明しなくても分かる、か。

 

 f:id:onsanai:20160819222324j:plain

 

サーマルグリスを除去したメインのチップ。カメラのレンズが映り込むくらい銀色に輝いている。

 

f:id:onsanai:20160819224547j:plain

 

 壊れていたのが電源回路のコイルの1個と思われる。ここがショートして電源の過電流保護が働いてパソコンが起動しなくなったようだ。保障は切れているのでELSAは悪くない。

 

f:id:onsanai:20160819222328j:plain

 

コイルが焼き切れたのがこれで2台目。このタイプがあと約10台あるので、続々と壊れるんだろうなあ(泣)。

 

ほかのいくつかの個体で、ヒートパイプからオイル漏れが確認できた。冷却性能がどこまで保たれているかは分からない。