ns-3でパケットのペイロードに対しデータを付加・取得する

NS-3でパケットに送信側でペイロードに独自データを付加して，受信側でそれを取得する方法についてのメモ．なおこの記事ではTCPソケットでパケットを送受信する想定です．

まず送信側で送るパケットのペイロードにメッセージを付加してパケットを生成・送信します．下の例ではstring文字列をuint8_tのポインタで保持してパケットのペイロードに付加しています．

std::ostringstream msg;
msg << "Hello, world!";
uint16_t packetSize = msg.str().length();
Ptr<Packet> packet = Create<Packet>((uint8_t *)msg.str().c_str(), packetSize);
m_socket->Send(packet);

受信側ではIPパケットの受信イベントを監視しておき，イベントが発生したら受信したパケットのペイロードからデータを取得します．

static void Rx(std::string context, Ptr<const Packet> packet, Ptr<Ipv4> ipv4,
               uint32_t interface) {
  Ptr<Packet> p = packet->Copy();

  Ipv4Header ipheader;
  TcpHeader tcpheader;
  p->RemoveHeader(ipheader);
  p->RemoveHeader(tcpheader);

  uint8_t *buffer = new uint8_t[p->GetSize()];
  p->CopyData(buffer, p->GetSize());
  std::string s = std::string(buffer, buffer + p->GetSize());
  std::cout << "Received:" << s << std::endl;
  delete buffer;
}

この例では

$ns3::Ipv4L3Protocol/Rx

をトレースしています．

ここでの注意は，受信側データを読むためには，送信したレイヤと同じレイヤまでパケットのヘッダを受信側で取り除く必要がある点です．今回はデータを付加したパケットをL4（TCP）レイヤで送信し，L3（IP）レイヤでパケットをトレースしているので，パケットからデータを読むにはIPヘッダとTCPヘッダを順に取り除きます．ですのでもしL2（MAC）レイヤでトレースしたのなら，まず該当するデータリンク層のヘッダから取り除く必要があります．

ns-3で指定したノードのIPv4アドレスを文字列で得る

ns-3でモデルを書いていると，今着目しているノードのアドレスが何なのかわからなくなるときがよくある．そんなときのためノードとインターフェースを指定してIPv4アドレスを文字列で得るためのメモ．

string getIpAddress (Ptr<Node> node)
{
  Ptr<Ipv4> ipv4 = node->GetObject<Ipv4> ();

  uint32_t ipvif = 1; // 0はループバックアドレス
  Ipv4InterfaceAddress ifAddr = ipv4->GetAddress (ipvif, 0);
  Ipv4Address addr = ifAddr.GetLocal ();

  uint32_t rawAddr = addr.Get ();
  ostringstream oss;
  oss << ((rawAddr >> 24) & 0xff) << "." << ((rawAddr >> 16) & 0xff) << "."
      << ((rawAddr >> 8) & 0xff) << "." << ((rawAddr >> 0) & 0xff);

  return oss.str ();
}

ノードのポインタを渡して，そのノードにインストールされているインターフェースの番号（0番がループバックアドレスになるので，今回はその次の番号．）を指定し，そのインターフェースに紐付いている最初のアドレスをuint32_tで取得して文字列に整形しています．

次のスクリプトで2台のノードをPointToPointネットワークで接続して，10.1.1.0/24を自動的に割り当てたあと，割り当てたアドレスを参照してみます．

int main (int argc, char **argv){
  const int NUM = 2;
  NodeContainer nodes;
  nodes.Create (NUM);

  PointToPointHelper pointToPoint;
  pointToPoint.SetDeviceAttribute ("DataRate", StringValue ("5Mbps"));
  pointToPoint.SetChannelAttribute ("Delay", StringValue ("2ms"));

  NetDeviceContainer devices;
  devices = pointToPoint.Install (nodes);

  InternetStackHelper stack;
  stack.Install (nodes);

  Ipv4AddressHelper address;
  address.SetBase ("10.1.1.0", "255.255.255.0");

  Ipv4InterfaceContainer interfaces = address.Assign (devices);

  for (int i = 0; i < NUM; i++){
    cout << getIpAddress (nodes.Get (i)) << endl;
  }
  return 0;
}

手元で実行すると以下のように出力されたのでとりあえずあってそうです．

10.1.1.1
10.1.1.2

ns-3で指定したパケットをロスさせる

ns-3で実験をするときにパケットを確率的にロスさせることができますが、ロスさせたいパケットを指定して落とすこともできるとわかったのでメモしておきます。

ロスするパケットを指定するにはエラーモデルのうちReceiveListErrorModelもしくはListErrorModelを使用します。前者は受信側で指定した順番のパケットを落とすのに対して、後者はパケットのUidを指定してロストさせるモデルです。前者との違いはUidを指定するため、受信側でパケットの到着順序が前後しても落とせるという点にあるようです。今回は前者のやりかたのメモです。

ロスの指定の方法は簡単で、ReceiveListErrorModelのインスタンスに対し、ロスさせたいパケット番号が入ったリストを設定します。そうしたらそのインスタンスを設定したいNetDeviceContainerの受信側ノードの"ReceiveErrorModel"Attributeに指定するだけです。

//ロスパケットの指定
//受信側で10,100,200,300番目に到着したパケットをロストさせる
std::list<uint32_t> packetErrorList = {10, 100, 200, 300};
Ptr<ReceiveListErrorModel> em = CreateObject<ReceiveListErrorModel> ();
em->SetList (packetErrorList);
//受信側ノードにエラーモデルを設定
devices.Get (1)->SetAttribute ("ReceiveErrorModel", PointerValue (em));

2ノード間をPointToPointネットワークで接続するシンプルなモデルで、受信側でパケットをロストさせる例を以下に示します。

#include "ns3/core-module.h"
#include "ns3/network-module.h"
#include "ns3/internet-module.h"
#include "ns3/point-to-point-module.h"
#include "ns3/applications-module.h"
#include <list>

#define PORT 10

using namespace ns3;

NS_LOG_COMPONENT_DEFINE ("ReceiveListErrorScript");

int num_pac = 0;

void inline MacRx (Ptr<const Packet> p){
  num_pac++;
};

void inline PhyRxDrop (Ptr<const Packet> p){
  std::cout << "Drop"
            << " " << num_pac++ << std::endl;
};

int main (int argc, char *argv[]){

  Time::SetResolution (Time::NS);
  LogComponentEnable ("UdpEchoClientApplication", LOG_LEVEL_ERROR);
  LogComponentEnable ("UdpEchoServerApplication", LOG_LEVEL_ERROR);

  NodeContainer nodes;
  nodes.Create (2);

  PointToPointHelper pointToPoint;
  pointToPoint.SetDeviceAttribute ("DataRate", StringValue ("10Mbps"));
  pointToPoint.SetChannelAttribute ("Delay", StringValue ("10ms"));
  
  NetDeviceContainer devices;
  devices = pointToPoint.Install (nodes);
  
  //ロスパケットの指定
  std::list<uint32_t> packetErrorList = {10, 100, 200, 300};
  Ptr<ReceiveListErrorModel> em = CreateObject<ReceiveListErrorModel> ();
  em->SetList (packetErrorList);
  devices.Get (1)->SetAttribute ("ReceiveErrorModel", PointerValue (em));
  
  InternetStackHelper stack;
  stack.Install (nodes);
  
  Ipv4AddressHelper address;
  address.SetBase ("10.1.1.0", "255.255.255.0");
  Ipv4InterfaceContainer interfaces = address.Assign (devices);
  
  UdpEchoServerHelper echoServer (PORT);
  ApplicationContainer serverApps = echoServer.Install (nodes.Get (1));
  serverApps.Start (Seconds (1.0));
  serverApps.Stop (Seconds (10.0));
  
  UdpEchoClientHelper echoClient (interfaces.GetAddress (1), PORT);
  echoClient.SetAttribute ("MaxPackets", UintegerValue (1000));
  echoClient.SetAttribute ("Interval", TimeValue (MilliSeconds (10)));
  echoClient.SetAttribute ("PacketSize", UintegerValue (1024));
  
  ApplicationContainer clientApps = echoClient.Install (nodes.Get (0));
  clientApps.Start (Seconds (2.0));
  clientApps.Stop (Seconds (10.0));
  
  //パケット受信・ロスのトレース
  std::string prefix = "/NodeList/1";
  std::string _macrx = prefix + "/DeviceList/*/$ns3::PointToPointNetDevice/MacRx";
  std::string _phyrxdrop = prefix + "/DeviceList/*/$ns3::PointToPointNetDevice/PhyRxDrop";
  
  Config::ConnectWithoutContext (_macrx, MakeCallback (&MacRx));
  Config::ConnectWithoutContext (_phyrxdrop, MakeCallback (&PhyRxDrop));
  
  Simulator::Run ();
  Simulator::Destroy ();
  
  return 0;
}

実行すると指定したパケットのときに指定したコールバック関数PhyRxDropが呼ばれます。

一つのNetDeviceに対してロスモデルの併用は不可能のようなので、こういったケースを模擬したい場合にはロスパケットを指定するモードにして、併用したケースを想定したパケットロスのリストを与えるのが楽だと思います。

tcコマンドのパケットロスモデル

tcコマンドを使ってネットワークをエミュレートするとき，パケットロスを設定することができます．tcコマンドのnetemには通常4種類のパケットロスのモデルが実装されています．かんたんに調べたのでメモしておきます．

各パケットロスモデル

ベルヌーイモデル (Bernoulli model)

　これは簡単で，ひとつのパケットのロス確率をpとするベルヌーイ試行にもとづく無記憶性のモデル．

（簡易）ギルバートモデル (Simplified Gilbert model)

　これは2つの状態を持つモデルで，連続したパケットロス（バーストロス）をモデル化できる．
　このモデルは2つのパラメータからなる．一つは受信良好状態からバーストロス状態に遷移する確率p，もう一つは反対にバーストロス状態から受信良好状態に遷移する確率r．

ギルバートモデル (Gilbert model)

　これは2つの状態を持つモデルで，連続したパケロス＋バーストロス中でも一定の割合で受信できる（ベルヌーイ分布），という状況をモデリングできる．
　このモデルは3つのパラメータからなる．一つは受信良好状態からバーストロス状態に遷移する確率p，もう一つは反対にバーストロス状態から受信良好状態に遷移する確率r．もう一つはバーストロス中であっても受信できる確率h（この場合，受信良好状態に遷移はしない）．
　h=0のとき，2.（簡易）ギルバートモデルと等価になる．

ギルバート・エリオットモデル (Gilbert-Elliot model)

　上の2つのモデルからなんとなく想像できる通り，3に加えて，受信良好状態であっても一定の割合でパケットをロス（ベルヌーイ分布）するケースを加えたモデル．
　このモデルは4つのパラメータからなる．受信良好状態からバーストロス状態に遷移する確率p，反対にバーストロス状態から受信良好状態に遷移する確率r．バーストロス中であっても受信できる確率h（この場合，受信良好状態に遷移はしない）．受信良好状態であってもパケットロスする確率k（この場合，バーストロス状態に遷移はしない）
　k=0のとき，3. ギルバートモデルと等価で，k=0, h=0のとき，2.（簡易）ギルバートモデルと等価になる．

tcコマンド上での指定

tcコマンドでは以下のようにパラメータを指定することで上記モデルを設定する．

tc qdisc add dev eth0 root netem loss p% r% 1-h% 1-k%

r, 1-h, 1-k を指定しなければ一番単純なベルヌーイモデルになり，全部指定すればギルバート・エリオットモデルになる．

プロクラステス分析

点群のレジストレーションや姿勢推定の論文で出てくるプロクラステス分析を知らなかったので調べて実装しました．

プロクラステス分析とは

プロクラステス分析（procrustes analysis，プロクラステス解析とも）は点同士の対応がとれた2つの点群に対して，並進・回転・一様なスケーリングの変換のもとで，点群間の二乗誤差が最小になるように重ねあわせする処理です．この記事では特に，回転行列が直交行列であるようなケースの直交プロクラステス分析(Orthogonal Procrustes analysis)について扱います．

方法

1. 2つの点群の重心を原点に位置合わせする

点群の重心が原点になるよう，からそれぞれの平均を引きます． \begin{equation} {\bf X}_{o} = {\bf X} - \overline{\bf X} \qquad {\bf Y}_{o} = {\bf Y} - \overline{\bf Y} \end{equation} 2. 2つの点群のスケールを正規化する

\begin{equation} {\bf X}_{norm} = \dfrac{{\bf X_o}}{||{\bf X}_o||_F} \qquad {\bf Y}_{norm} = \dfrac{{\bf Y_o}}{||{\bf Y}_o||_F} \end{equation} 3. 回転に二乗誤差が最小になるように回転する

重心を原点に位置合わせし，スケールを正規化した2つの点群を改めてとします．すると，を回転して誤差が最小になるようにに重ねあわせる時の回転行列は以下の式で求めることができます． \begin{equation}
{\bf R} = argmin_{\bf \Omega}{|| {\bf \Omega A} - {\bf B} ||}_F \quad s.t. \quad {\bf \Omega}^{T} {\bf \Omega} = I \end{equation}

この式は
\begin{equation} ||{\bf \Omega A} - {\bf B}||^2_F = ||{\bf A}||^2_F + ||{\bf B}||^2_F - 2 \langle {\bf \Omega A}, {\bf B}\rangle \end{equation} と変形でき、これを最小化するには右辺第三項を最大化します．次にこの項を以下のように変形します．

\begin{align*} \langle{\bf \Omega A}, {\bf B}\rangle &= tr({\bf A^T \Omega^T B}) = tr({\bf B^T \Omega A}) = tr({\bf \Omega A B^T}) \newline &= tr({\bf \Omega USV^T})= tr({\bf SV^T\Omega U})
\end{align*} となります． （ただしと特異値分解した．）

と置くと，この右辺は直交行列の積であることからは直交行列になり，式は以下のように変形できます．
\begin{equation} tr({\bf SV^T\Omega U}) = tr({\bf SX}) = \sum_i s_{ii}x_{ii}
\end{equation} は特異値であるため非負になるので，この式を最大化する直交行列は明らかにの時です． よってであるため の時，当初の式が最小化されるのでこれが求めたい回転行列になります．

実装

実装は↓にあるけど，ためしに自分でも実装してみました． stackoverflow.com

自分の実装

import numpy as np

def procrustes(x, y):
    """
    Orthogonal procrustes analysis

    Inputs:
    ------------
    X, Y    
        matrices of target and input coordinates. 

    Outputs
    ------------
    Z
        the matrix of transformed X-values
    """
    mux = np.mean(x, axis=0)
    muy = np.mean(y, axis=0)
    X = x - mux
    Y = y - muy
    X_norm = np.linalg.norm(X) # frobenius norm
    Y_norm = np.linalg.norm(Y) # frobenius norm
    normalized_X = X / X_norm
    normalized_Y = Y / Y_norm
    
    u, w, vt = np.linalg.svd(np.dot(normalized_X, normalized_Y.T))
    
    R = np.dot(vt.T, u.T)
    scale = w.sum() * X_norm / Y_norm
    Z = Y_norm * w.sum() * np.dot(R, normalized_X) + muy 
    
    return Z

二次元の点群に対して回転・並進・スケーリングした上で各点にガウスノイズを加えてできた点群と，もとの点群でプロクラステス分析した結果が次の図です．

補足

プロクラステス分析は2点群間の対応が取れているのに対して，ICP(Iterative Closest Point)は一般に揃えたい点群間の対応が取れていない場合や点数が異なるケースを扱っている（大規模な点群の一部分に位置合わせしたい，など）という違いがあるようです．

参考資料

[1] Ross, A. "Procrustes analysis." Course report, Department of Computer Science and Engineering, University of South Carolina (2004).
[2] Akça, M. Mehmet, D. "Generalized procrustes analysis and its applications in photogrammetry." (2003).

NNablaでVGG16を実装

ソニー社内で使用されているというDNNフレームワークがNNablaという名前でOSS化されたので使ってみました．練習として，画像分類で代表的なモデルのVGG16-netを実装してみました．

import nnabla as nn

import nnabla.functions as F
import nnabla.parametric_functions as PF
import nnabla.solvers as S

import numpy as np
from collections import OrderedDict
import cv2

class VGG16(object):
    def __init__(self, input_tensor, weights='imagenet'):
        self.batch_size = input_tensor.shape[0]
        self.layers = OrderedDict()
        
        # Input
        self.x = self.layers['input'] = input_tensor
        
        # Convolution layers
        h = self.layers['block1_conv1'] = PF.convolution(self.x, 64, (3, 3), pad=(1, 1), stride=(1, 1), name='block1_conv1')
        h = F.relu(h)
        h = self.layers['block1_conv2'] = PF.convolution(h, 64, (3, 3), pad=(1, 1), stride=(1, 1), name='block1_conv2')
        h = F.relu(h)
        h = self.layers['block1_pool1'] = F.max_pooling(h, (2, 2), stride=(2, 2))
        
        h = self.layers['block2_conv1'] = PF.convolution(h, 128, (3, 3), pad=(1, 1), stride=(1, 1), name='block2_conv1')
        h = F.relu(h)
        h = self.layers['block2_conv2'] = PF.convolution(h, 128, (3, 3), pad=(1, 1), stride=(1, 1), name='block2_conv2')
        h = F.relu(h)
        h = self.layers['block2_pool1'] = F.max_pooling(h, (2, 2), stride=(2, 2))
        
        h = self.layers['block3_conv1'] = PF.convolution(h, 256, (3, 3), pad=(1, 1), stride=(1, 1), name='block3_conv1')
        h = F.relu(h)
        h = self.layers['block3_conv2'] = PF.convolution(h, 256, (3, 3), pad=(1, 1), stride=(1, 1), name='block3_conv2')
        h = F.relu(h)
        h = self.layers['block3_conv3'] = PF.convolution(h, 256, (3, 3), pad=(1, 1), stride=(1, 1), name='block3_conv3')
        h = F.relu(h)
        h = self.layers['block3_pool1'] = F.max_pooling(h, (2, 2), stride=(2, 2))
        
        h = self.layers['block4_conv1'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block4_conv1')
        h = F.relu(h)
        h = self.layers['block4_conv2'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block4_conv2')
        h = F.relu(h)
        h = self.layers['block4_conv3'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block4_conv3')
        h = F.relu(h)
        h = self.layers['block4_pool1'] = F.max_pooling(h, (2, 2), stride=(2, 2))
        
        h = self.layers['block5_conv1'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block5_conv1')
        h = F.relu(h)
        h = self.layers['block5_conv2'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block5_conv2')
        h = F.relu(h)
        h = self.layers['block5_conv3'] = PF.convolution(h, 512, (3, 3), pad=(1, 1), stride=(1, 1), name='block5_conv3')
        h = F.relu(h)
        h = self.layers['block5_pool1'] = F.max_pooling(h, (2, 2), stride=(2, 2))
        
        # Flatten
        h = F.transpose(h, (0, 2, 3, 1)) # (batch, channels, h, w) -> (batch, h, w, channels)
        h = self.layers['flatten'] = F.reshape(h, (self.batch_size, np.prod(h.shape[1:])))
        
        # Fully-Connected layers
        h = self.layers['fc1'] = PF.affine(h, 4096, name='fc1')
        h = F.relu(h)
        h = self.layers['fc2'] = PF.affine(h, 4096, name='fc2')
        h = F.relu(h)
        h = self.layers['fc3'] = PF.affine(h, 1000, name='fc3')
        
        # Prediction
        output = self.layers['prob'] = F.softmax(h)
        
        self.output = output
    
        
    def summary(self):
        layer_len = len(self.layers.keys())
        maxlen = max(map(len, self.layers.keys()))
        print("{} layers.".format(layer_len))
        for k in self.layers:
            print("{} layer: {}".format(k.ljust(maxlen), str(self.layers[k].shape)))


x = nn.Variable([1, 3, 224, 224])
vgg16 = VGG16(x)

ここから学習するには学習データを用意したうえで，lossとSolverを定義してforwardしてbackwardしてという手順になるのですが，あいにく肝心のimagenetのデータを持っていません．なので今回はKerasで使用されている学習済みの重みを移植することにしました．
まず以下のようにKeras側でpre-trainedなVGG16-netを構築し，学習済みのweightをpickleで吐き出します。

import pickle

model=  VGG16(weights="imagenet")
weights = model.get_weights()

mat = []
for l in weights:
    if len(l.shape) is 4:
        l = l.transpose()
    mat.append(l)

with open('keras_nn.pickle', 'wb') as f:
    pickle.dump(mat, f)

pickle化したweightを，以下のコードでnnablaのVGG16-netに渡しています．（ここはもっといいやり方があると思う）

params = nn.get_parameters()
with open('keras_nn.pickle', 'rb') as f:
    mats = pickle.load(f)
    
idx = 0
suffix = ['/conv/W', '/conv/b', '/affine/W', '/affine/b']
for k in alex.layers:
    for suf in suffix:
        w_name = k + suf
        if w_name not in params:
            continue
        if suf is '/conv/W':
            mats[idx] = mats[idx].transpose((0,1,3,2))
        params[w_name].d = np.float32(mats[idx])
        idx += 1

ここまででnnablaでVGG16-netの構築と学習済みの重みを読み込みました．次に実際に画像を読み込んでクラス分類してみます．Kerasのモデルは色チャンネルをBGRの順で読みこんでいるのでopencvのimreadで読み込むと良いです．

def prepare_img(img):
    img = img.astype(np.float32)
    img = cv2.resize(img, (224, 224))
    img[:, :, 0] -= 103.939
    img[:, :, 1] -= 116.779
    img[:, :, 2] -= 123.68
    img = np.expand_dims(img, axis=0)

    img = img.transpose((0, 3, 2, 1))

    return img

img = cv2.imread('input.jpg', 1)
preproc_img = prepare_img(img)
x.d = preproc_img
vgg16.output.forward()
label = np.argsort(vgg16.output.d[0])[::-1][:5]
for i in label:
    print(i, vgg16.output.d[0][i])

出力は以下のようになりました．

792 0.985962 //shovel
587 0.00602465 //hammer
813 0.00355319 //spatula
596 0.00200536 //hatchet
462 0.000546574 //broom

847 0.982786 //tank
471 0.0112128 //cannon
744 0.00218588 //projectile, missile
408 0.00189686 //amphibian
657 0.00173452 //missile

とりあえずうまく行っているかな?