どすえのブログ

京都在住プログラマーの開発ブログ。バイクとキャンプが趣味。

Fortranによる3次元線形補間 | 3D Linear Interpolation

Fortranで3次元の線形補完ルーチンを書く機会があったのでメモ。1次元のケースから確認し、2次元、3次元に拡張する流れで進める。

1次元の線形補間

1次元グリッド x_1, x_2, ... x_n上で値 v_1, v_2, ... v_nがそれぞれ定義されているとする.

ここでは,  x_1 x_2の間にある x_iにおける値 v_iを求めてみる.

まず距離として以下の量を定義する.

 
\begin{align}
l_1 &= x_i - x_1 \\
l_2 &= x_2 - x_i \\
l &= l_1 + l_2
\end{align}

 (x_i,v_i) (x_1,v_1) (x_2,v_2)を結んだ直線上に位置しているので,  v_i v_1 v_2の距離に逆比例する重み付き平均で求められる. 式で書けば

 
\begin{align}
v_i &= v_1 \frac{l_2}{l} + v_2 \frac{l_1}{l} \\
&= v_1 \frac{x_2-x_i}{x_2 - x_1} + v_2 \frac{x_i - x_1}{x_2-x_1}
\end{align}

となる. このイメージを持ったまま多次元に拡張すると分かりやすい.

f:id:dosuex:20210526182936p:plain
線形補間のイメージ図

2次元の線形補間

2次元に拡張した場合を考える.
まずx軸方向で1次元線形補間を行い(図中の灰色点), 次にy軸方向で1次元線形補間を行う(図中の赤色点)と考えれば良い. 式で書くと

 
\begin{align}
v_{ii} &= \frac{y_2 - y_i}{y_2 - y_1}\left( v_{11} \frac{x_2-x_i}{x_2 - x_1} + v_{21} \frac{x_i - x_1}{x_2-x_1} \right) + \frac{y_i - y_1}{y_2 - y_1}\left( v_{12} \frac{x_2-x_i}{x_2 - x_1} + v_{22} \frac{x_i - x_1}{x_2-x_1} \right) \\
&= \frac{1}{S} \{ (x_2-x_i)( y_2 - y_i )v_{11} + (x_i - x_1)( y_2 - y_i )v_{21} \\
&\quad + (x_2-x_i)(y_i - y_1)v_{12} + (x_i - x_1)(y_i - y_1)v_{22} \}
\end{align}

ただし S=(x_2-x_1)(y_2-y_1)としている. 式を眺めると添字の関係性に規則性があることがわかる.

f:id:dosuex:20210526192038p:plain
2次元線形補間イメージ図

3次元の線形補間

3次元の図を書くのが大変なのであとは規則性に従って公式を導く.

 
\begin{align}
v_{ii} = \frac{1}{V} \{ (x_2-x_i)(y_2 - y_i )(z_2 - z_i )v_{111} + (x_2-x_i)(y_i - y_1 )(z_2 - z_i )v_{121}\\ + (x_i-x_1)(y_i-y_1)(z_2-z_i)v_{221}
+ (x_i-x_1)(y_2-y_i)(z_2-z_i)v_{211}\\ + (x_i-x_1)(y_2-y_i)(z_i-z_1)v_{212} + (x_i-x_1)(y_i-y_1)(z_i-z_1)v_{222}\\
+ (x_2-x_i)(y_i-y_1)(z_i-z_1)v_{122} + (x_2-x_i)(y_2-y_i)(z_i-z_1)v_{112}
 \}
\end{align}

ただし V=(x_2-x_1)(y_2-y_1)(z_2-z_1)としている.

コード

Fortranによる3次元の線形補間アルゴリズムを載せます.

module interpol

    implicit none
    
    type InterpolCoeffs3D
        integer :: i = 0
        integer :: j = 0
        integer :: k = 0
        real(8) :: b111 = 0.d0
        real(8) :: b121 = 0.d0
        real(8) :: b112 = 0.d0
        real(8) :: b122 = 0.d0
        real(8) :: b211 = 0.d0
        real(8) :: b212 = 0.d0
        real(8) :: b221 = 0.d0
        real(8) :: b222 = 0.d0
    end type InterpolCoeffs3D
    
    contains
    
    subroutine xyz_interpol3D_coeff(xmin,dx,nx,ymin,dy,ny,zmin,dz,nz,xout,yout,zout,c,err)
        real(8), intent(in)           :: xmin
        real(8), intent(in)           :: dx
        integer, intent(in)           :: nx
        real(8), intent(in)           :: ymin
        real(8), intent(in)           :: dy
        integer, intent(in)           :: ny
        real(8), intent(in)           :: zmin
        real(8), intent(in)           :: dz
        integer, intent(in)           :: nz
        real(8), intent(in)           :: xout
        real(8), intent(in)           :: yout
        real(8), intent(in)           :: zout
        type(InterpolCoeffs3D), intent(out) :: c
        integer, intent(out), optional      :: err
    
        real(8) :: x1, x2, y1, y2, z1, z2, xp, yp, zp, dV
        integer :: i, j, k, err_status
    
        err_status = 1
        xp = max(xout,xmin)
        yp = max(yout,ymin)
        zp = max(zout,zmin)
        i = floor((xp-xmin)/dx)+1
        j = floor((yp-ymin)/dy)+1
        k = floor((zp-zmin)/dz)+1
    
        if ((((i.gt.0).and.(i.le.(nx-1))).and.((j.gt.0).and.(j.le.(ny-1)))).and.((k.gt.0).and.(k.le.(nz-1)))) then
            x1 = xmin + (i-1)*dx
            x2 = x1 + dx
            y1 = ymin + (j-1)*dy
            y2 = y1 + dy
            z1 = zmin + (k-1)*dz
            z2 = z1 + dz
            dV = dx*dy*dz
    
            c%b111 = ((x2 - xp) * (y2 - yp) * (z2 - zp))/dV
            c%b121 = ((x2 - xp) * (yp - y1) * (z2 - zp))/dV
            c%b221 = ((xp - x1) * (yp - y1) * (z2 - zp))/dV
            c%b211 = ((xp - x1) * (y2 - yp) * (z2 - zp))/dV
            c%b212 = ((xp - x1) * (y2 - yp) * (zp - z1))/dV
            c%b222 = ((xp - x1) * (yp - y1) * (zp - z1))/dV
            c%b122 = ((x2 - xp) * (yp - y1) * (zp - z1))/dV
            c%b112 = ((x2 - xp) * (y2 - yp) * (zp - z1))/dV
            c%i = i
            c%j = j
            c%k = k
            err_status = 0
        endif
    
        if(present(err)) err = err_status
    
    end subroutine xyz_interpol3D_coeff
    
    end module interpol
    
    program test
        use interpol
        implicit none

        real(8) :: xmin,dx,ymin,dy,zmin,dz,xout,yout,zout,rnd,vout
        integer :: nx,ny,nz,err,ix,iy,iz,i,j,k
        type(InterpolCoeffs3D) :: c
        real(8),allocatable :: v(:,:,:)

        xmin = 1.0d0
        dx   = 0.5d0
        ymin = 2.0d0
        dy   = 0.4d0
        zmin = 3.0d0
        dz   = 0.6d0
        nx   = 10
        ny   = 15
        nz   = 5

        allocate(v(nx,ny,nz))

        do ix = 1,nx
            do iy = 1,ny
                do iz = 1,nz
                    call random_number(rnd)
                    v(ix,iy,iz) = rnd
                end do
            end do
        end do

        call xyz_interpol3D_coeff(xmin,dx,nx,ymin,dy,ny,zmin,dz,nz,xout,yout,zout,c,err)

        i = c%i
        j = c%j
        k = c%k
        vout = c%b111*v(i,j,k) + c%b121*v(i,j+1,k) + c%b221*v(i+1,j+1,k) + c%b211*v(i+1,j,k) + &
               c%b212*v(i+1,j,k+1) + c%b222*v(i+1,j+1,k+1) + c%b122*v(i,j+1,k+1) + c%b112*v(i,j,k+1)

        write(*,*) vout

    end program test
    

非線形カルマンフィルタ (2) - アンサンブルカルマンフィルタ -

目次

本シリーズについて

非線形現象に対するカルマンフィルタを本で勉強したのでメモも兼ねてまとめます.

前回の記事 dosuex.hatenablog.com では状態遷移、観測がともに線形である(行列演算で記述できる)ケースに適用できるカルマンフィルタを紹介しました。

今回は非線形システムを扱うことができるアンサンブルカルマンフィルタ(Ensemble Kalman Filter; EnKF)を紹介します. アンサンブルカルマンフィルタは, モンテカルロ近似を用いて, 非線形のシステムモデルをそのまま使って分散共分散行列の推定を行います。

コード( Jupyter Notebook )はこちらです.

github.com

今回はシステムモデルに関しては非線形モデルを, 観測モデルに関しては線形モデルを仮定する.

 
\begin{align}
x_{t} &= f_t(x_{t-1}, v_t) \tag{1} \\
y_t &= H_t x_t + w_t \tag{2} \\
\end{align}

x_t状態ベクトル,  y_t は観測データを表す. 観測ノイズは w_t \sim N(0,R_t)としてガウス分布からサンプルされるノイズを仮定する. 一方, システムノイズ v_t \sim q_t(v_t)は任意の分布で良い.

EnKFではアンサンブル近似を用いて分布の近似を行う. アンサンブル近似とは, 以下の図のようにある確率変数 xの確率分布を N個のサンプルの集合を用いて近似することである.

f:id:dosuex:20210517222636p:plain

このサンプルの集合 \{x^{(i)}\}_{i=1}^N をアンサンブルを呼ぶ. またそれぞれの x^{(i)} のことを粒子と呼ぶ.

EnKFでは, (1)式, (2)式で与えられる状態空間モデルのもとで, カルマンフィルタと同様に1期先予測とフィルタの手順を繰り返す.

カルマンフィルタでは平均値と共分散行列を逐次的に計算したが, EnKFでは予測分布とフィルタ分布のアンサンブルを逐次的に計算することになる.

これを改めて式で表す.

時刻t-1でのフィルタ分布p(x_{t-1}|y_{1:t-1})を近似するN個の粒子から構成されるアンサンブル \{x^{(i)}_{t-1|t-1}\}_{i=1}^N が準備されているものとする. これは

 
\begin{align}
p(x_{t-1}|y_{1:t-1}) \approx \frac{1}{N}\sum_{i=1}^N \delta (x_{t-1} - x^{(i)}_{t-1|t-1})
\end{align}

と書ける. このフィルタアンサンブルを使って時刻tでの観測y_tを取り込んで

 
\begin{align}
p(x_{t}|y_{1:t}) \approx \frac{1}{N}\sum_{i=1}^N \delta (x_{t} - x^{(i)}_{t|t})
\end{align}

となる \{x^{(i)}_{t|t}\}_{i=1}^N を求める一連の手続きが, EnKFである.

以下では, 添字の  t|t-1 は時刻  t-1 の情報を用いて時刻tの状態を推定する操作(予測)を表し,  t|tは現在の時刻の情報を用いて状態を最適化する操作(フィルタリング)を表す.

またプライム記号\primeは行列の転置を表す.

アンサンブルカルマンフィルタ

  1. 初期状態のアンサンブル \{x_{0|0}^{(i)}\}_{i=1}^N を生成する. t=1とする.

  2. 以下を行う.

    1. 1期先予測

      • システムノイズのアンサンブル \{ v_t^{(i)} \}_{i=1}^N を生成する.

      • iに対して x_{t|t-1}^{(i)} = f_t (x_{t-1|t-1}^{(i)}, v_t^{(i)}) を計算して, 1期先予測アンサンブルを得る.

    2. フィルタ

      • システムノイズのアンサンブル \{ w_t^{(i)} \}_{i=1}^N

      • 分散共分散行列 \hat{V}_{t|t-1} ,  \hat{R}_t ならびに \hat{K}_t を計算する

      • iに対して,  x_{t|t}^{(i)} = x_{t|t-1}^{(i)} + \hat{K}_t (y_t + \check{w}_t^{(i)} - H_t x_{t|t-1}^{(i)}) を計算し, フィルタアンサンブル \{x_{t|t}^{(i)}\}_{i=1}^N を得る.

3: t=Tならば停止. それ以外は t\leftarrow t+1として2に戻る.

以上の手続きで使用した変数の定義を示す.

  • 状態変数ベクトルの偏差
 
\begin{align}
\check{x}_{t|t-1}^{(i)} = x_{t|t-1}^{(i)} - \frac{1}{N} \sum_{j=1}^N x_{t|t-1}^{(j)}
\end{align}
  • 状態変数ベクトルの分散共分散行列
 
\begin{align}
\hat{V}_{t|t-1} = \frac{1}{N-1} \sum_{j=1}^N x_{t|t-1}^{(j)} x_{t|t-1}^{(j)^{\prime}}
\end{align}
  • 観測ノイズの偏差
 
\begin{align}
\check{w}_{t}^{(i)} = w_{t}^{(i)} - \frac{1}{N} \sum_{j=1}^N w_{t}^{(j)}
\end{align}
  • 観測ノイズの分散共分散行列
 
\begin{align}
\hat{R}_{t} = \frac{1}{N-1} \sum_{j=1}^N w_{t}^{(j)} w_{t}^{(j)^{\prime}}
\end{align}
  • カルマンゲイン
 
\begin{align}
\hat{K}_{t} = \hat{V}_{t|t-1} H_t^{\prime} (H_t \hat{V}_{t|t-1} H_t^{\prime} + \hat{R}_t)^{-1}
\end{align}

適用

非線形の減衰を受けた振動子の運動を記述するVan der Pol 方程式

 
\begin{align}
\frac{d^{2} x}{d t^{2}}-\mu\left(1-x^{2}\right) \frac{d x}{d t}+x=0
\end{align}

にEnKFを適用する. 減衰の強さを制御するパラメータである \muを未知数とし, このパラメータの真値を同定することで, 真の力学モデルを同定する問題を解く.

下図に適用した結果を示す. 赤線が真値であり, 緑色の三角が観測ノイズの乗った観測である. EnKFによる推定結果である青線が、次第に真値に近づいていく様子が見られる. f:id:dosuex:20210518015005p:plain

こちらは減衰パラメータの時間ステップごとの推定結果を示している. 青線がそれぞれのアンサンブル粒子(状態変数ベクトル)の保持する減衰パラメータに相当する. 全ての粒子が次第に正しい減衰パラメータの値に収束していく様子が見られる. f:id:dosuex:20210518224912p:plain

参考書籍

非線形カルマンフィルタ (1) - 線形カルマンフィルタ -

目次

本シリーズについて

非線形現象に対するカルマンフィルタを本で勉強したのでメモも兼ねてまとめます.

観測データに基づいて, 線形確率システムの状態ベクトルを逐次的に最適化・推定するアルゴリズムとして, カルマンフィルタが有名です. カルマンフィルタは状態ベクトルの時間発展と, 状態ベクトルの観測に線形性を仮定しているため, 非線形の現象に適用するためには工夫が必要となります. 本シリーズでは線形カルマンフィルタから順を追って非線形カルマンフィルタの実装へ進んでいきます.

今回は状態遷移、観測がともに線形である(行列演算で記述できる)ケースを考えます.

TL;DR

コード( Jupyter notebook )はこちらです.

github.com

状態空間モデルで表される線形確率システム

 
\begin{align}
x_{t+1} &= F_t x_t + G_t w_t \\
y_t &= H_tx_t + v_t\\
\end{align}

を考えます。ここで,  x_t\in \mathbb{R}^n  状態ベクトル,  y_t \in \mathbb{R}^p  は観測ベクトル,  w_t \in \mathbb{R}^m  ,  v_t \in \mathbb{R}^p ガウスノイズベクトルです. また F_t \in \mathbb{R}^{n\times n} は遷移行列,  G_t \in \mathbb{R}^{n \times m} は駆動行列,  H_t \in \mathbb{R}^{p \times n} は観測行列と呼ばれます. さらに, x_t, w_t,  v_t それぞれに対する誤差共分散行列としてP,Q,Rを導入します.

取り急ぎ, 下式たちを実装すればカルマンフィルタができます. 添字の t / t-1 は時刻  t-1 の情報を用いて時刻tの状態を推定する操作(予測)を表し, t/tは現在の時刻の情報を用いて状態を最適化する操作(フィルタリング)を表します. また t/N は時刻Nまで予測が完了している状態での時刻tの最適化操作(スムージング)を表します.

カルマンフィルタ

1:初期値を  \hat{x}_{0/-1}=\bar{x}_0 , P_{0/-1}=P_0 とおき,  t=0 とする.
2:観測更新ステップ

 A:カルマンゲイン

 
K_t = P_{t/t-1}H_t^T [H_t P_{t/t-1} H_t^T + R_t]^{-1}

 B:フィルタ推定値

 \hat{x}_{t/t} = \hat{x}_{t/t-1} + K_t [y_t - H_t \hat{x}_{t/t-1}]

 C:フィルタ推定誤差共分散行列

 
P_{t/t} = P_{t/t-1} - K_t H_t P_{t/t-1}

3:時間更新ステップ

 A:1段先予測値

 
\hat{x}_{t+1/t} = F_t  \hat{x}_{t/t}

 B:予測誤差共分散行列

 
P_{t+1/t} = F_t P_{t/t} F_t^T + G_t Q_t G_t^T

4: t \leftarrow t+1 としてステップ2に戻る.

カルマンスムーザ


\begin{align}
\hat{x}_{t/N} &= \hat{x}_{t/t} + C_t(\hat{x}_{t+1/N} - \hat{x}_{t+1/t})\\
C_t &= P_{t/t} F_t^T P_{t+1/t}^{-1} \\
P_{t/N} &= P_{t/t} + C_t [P_{t+1/N}-P_{t+1/t}]C_t^T
\end{align}

実装

動的な線形システムの逐次予測をカルマンフィルタを使って行います. 参考資料のp.56にある例の通り, 人工衛星の回転運動を線形近似した4次元システムを考えます.


\begin{align}
\dot{x}_t^{(1)} &= x_t^{(2)}\\
\dot{x}_t^{(2)} &= x_t^{(3)} + x_t^{(4)}\\
\dot{x}_t^{(3)} &= 0\\
\dot{x}_t^{(4)} &= -0.5 x_t^{(4)} + \xi_t
\end{align}

ただし,  {x}_t^{(1)} : 人工衛星の姿勢角度,  {x}_t^{(2)} : 角速度,  {x}_t^{(3)} : 角加速度の平均値成分,  {x}_t^{(4)} : 角加速度のランダム成分であり,  \xi_t は平均0のガウスノイズです. サンプル間隔  \Delta =1.0 で離散化すると時間発展は以下で記述されます.


x_{t+1} = \left(
\begin{array}{ccc}
1 & 1 & 0.5 & 0.5\\
0 & 1 & 1 & 1\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 0.606
\end{array}
\right) x_t + \left(
\begin{array}{ccc}
0\\
0\\
0\\
1
\end{array}
\right) w_t

ただし, ガウスノイズ w_t, v_t の分散は Q=0.64\times 10^{-2} ,  R=1 とします. 状態ベクトルの初期値および推定値の初期値を


x_0 = \left(\begin{array}{ccc}
1.25\\
0.06\\
0.01\\
-0.003
\end{array}
\right),

\hat{x}_{0/-1} = \left(\begin{array}{ccc}
0\\
0\\
0\\
0
\end{array}
\right),

P_{0/-1} = \rm{diag}[10,10,10,10 ]

としてシミュレーションを行います. 衛星の姿勢角推定の結果を下図に示します. f:id:dosuex:20210516124654p:plain

緑線のフィルタ推定値は観測値の影響を大きく受けていますが, 赤線のスムージング推定値は姿勢角の真の動きに迫っている様子が見えます.

参考資料

サンフランシスコで後ろから殴られてヘッドホンを盗られた話

カメラロールを見返していたら、こんな事件があったなと思い出したので振り返りブログに。

 

アメリカ西海岸を単身旅行中にサンフランシスコで殴り飛ばされてヘッドホンを盗まれた話です。

時は2018年2月、2週間のアメリカ西海岸旅行に出発しました。サンフランシスコ、ヨセミテ国立公園、ロサンゼルス、サンノゼを回る計画でした。事件は到着してわずか3日目のサンフランシスコで起こりました。

3日目はカリフォルニア屈指の景勝地であるヨセミテ国立公園に行く予定でした。移動用の長距離バスを予約していたので、朝5時にドミトリーの宿を出発してバス乗り場に向かっていました。サンフランシスコは大手町のようなビル群が立ち並ぶビジネス街なので、早朝はひと気もなく暗いです。冬はサンフランシスコでも寒く、耳あてがわりにオーバーイヤーヘッドホンをつけて歩いていました。

 アメリカではwi-fi以外でスマホをインターネットに繋げないので、地球の歩き方についているマップを使ってバス停まで歩くことにしました。地球の歩き方には旅人への注意が書いてあり、サンフランシスコの欄は事前にチェック済みです。それによると、サンフランシスコは綺麗な街だが、中心部にゴロツキの溜まり場となっているエリアがあり、ひったくりの被害は意外に多いので注意と書かれていました。

「こんな綺麗な街でそんなことあるのか?」と思いながらも、暗い路地(路地といっても車道二車線と歩道がある通り)にふとさしかかりました。今までとは明らかに雰囲気が違う場所だったので、ひったくりとかはこういう場所で起こったりするのかな。とか思いスマホで周囲の写真を撮りながら歩き続けました。

f:id:dosuex:20210515225147p:plain

その時撮った写真。なぜライドシェアを撮っているかは不明。

一通り撮り終えてスマホをポケットにしまったまさにその時。後ろから「タッタッタッタ」と走ってくる音が聞こえました。状況的にどう考えてもやばいのですが、ここは生粋の日本人である私、平和ボケは世界トップレベルです。

「朝ランかな?やっぱサンフランシスコ民は意識高いな。」

とか呑気なことを思いながら、後ろを振り返ることなく歩き続けました。次の瞬間、左顎に強烈な衝撃を感じ、そのまま地面に転がり倒れました。

何が起きたのかさっぱりわかりませんでしたが、後ろを振り返るとムキムキでゴリッゴリの黒人が私に向かって大声で吠えていました。右手には私がつけていたヘッドホンが握られていました。彼はそれが欲しかったのだろうと思いヘッドホン以外の地面に散乱した荷物を抱えて逃げました。

本当に急なアクシデントが起きると、人間は痛みとか焦りとか恐怖とかは感じないんですね。その後数分間はやけに冷静というか、ちょっと世界がスローに感じました。これがアドレナリンでしょうか。

バス亭についたくらいで左頬と目の横あたりが強烈に痛み出し、頬はパンパンに腫れているのがはっきりとわかりました。

 

それでもここは海外、泣きつく相手もいないのでそのまま旅程を続行しヨセミテへ向かいました。人生でも最大級に孤独感と無力感を感じた1日でした。

f:id:dosuex:20210515225750p:plain

ただしヨセミテは本当に素晴らしい場所でこんなことがあった後でも最高でした。また行きたい。

夜、ヨセミテの宿について鏡を見ると頬が腫れすぎて直角になっており、さらにアザで黄色くなっていて見た目はさしずめ食パンマン。

殴られただけでナイフで刺されたりしなくてよかった、盗られたのがヘッドホンだけでよかった、とか思いながら泥のように眠りました。

 

いやーそれにしても、ガリガリのアジア人が地球の歩き方持って、耳をヘッドホンで塞いで暗い道を一人で歩くってカモすぎましたね。しかも走ってきても振り返らなかった。笑

地球の歩き方は表紙が特徴的で多くの日本人観光客が持っているので、堂々と持ち歩いてるとひったくりのターゲットにされやすいみたいです。みなさんも気をつけてください。

threadingとQueueを使ったDjangoでのマルチスレッド処理

Djangoで、リクエストが来たら発火する処理が時間のかかるものだった場合、処理自体は裏で走らせておいて、とりあえずレスポンスを返しておきたいというケースがあると思います。

例えば、ユーザー登録が完了した際に確認メールを送信するといった作業です。メール送信に1秒かかるアプリケーションに、秒間1000アクセスきた場合

  • メール送信のネットワークI/Oの間にアプリケーションがロックされない。
  • メール送信のジョブを裏で走らせて、アクセスがあった順序で実行していく。

ということが求められるでしょう。

Djangoでメール送信機能を作る機会があったので、メモ。

基本動作の確認

まずはqueueのリファレンスにあるサンプルコードを少し改変してqueueとthreadingを使った処理を確認。

import threading
import queue
import time

q = queue.Queue()

def worker():
    while True:
        item = q.get()
        print(f'Working on {item}')
        time.sleep(1)
        print(f'Finished {item}')
        q.task_done()

# turn-on the worker thread
threading.Thread(target=worker, daemon=True).start()

# send thirty task requests to the worker
for item in range(5):
    print("Put item", item)
    q.put(item)
print('All task requests sent\n', end='')

# block until all tasks are done
q.join()
print('All work completed')

最初にQueueインスタンスとしてqグローバル変数で宣言します。このキューに実行すべきジョブが溜まります。次にジョブの実行主体であるworkerメソッドを定義します。workerではジョブを取得してprintし、1秒待つという処理を行います。while Trueで常に待機状態にありますが、q.get()が返り値を戻した時だけその下の処理が実行されます。 次にthreading.Thread()で新規スレッドを立ち上げます。targetに実行したいメソッドを渡します。 あとはキューにジョブとしてitemを追加するとworkerが走ります。q.join()qからgetされた全てのジョブに対して完了宣言q.task_done()が行われるまでそれ以降の処理を止める役割をします。

出力

Put item 0
Put item 1
Put item 2
Put item 3
Put item 4
All task requests sent
Working on 0
Finished 0
Working on 1
Finished 1
Working on 2
Finished 2
Working on 3
Finished 3
Working on 4
Finished 4
All work completed

出力を見ると最初にジョブの追加が全て完了し、追加した順にジョブが実行されており、非同期的にジョブ管理・実行が行われていることがわかります。この要領でメール送信も実装できそうです。試してみます。

Djangoでのタスクの非同期処理

まず、テスト用のDjango APIを作成します。

$ django-admin startproject mysite
$ cd mysite/
$ python manage.py startapp api

そして以下のようにコードを書き加えてください。

mysite/settings.py

(中略)
INSTALLED_APPS = [
    'django.contrib.admin',
    'django.contrib.auth',
    'django.contrib.contenttypes',
    'django.contrib.sessions',
    'django.contrib.messages',
    'django.contrib.staticfiles',
    'api',
    'rest_framework'
]
(中略)

mysite/urls.py

from django.contrib import admin
from django.urls import path
from api import views

urlpatterns = [
    path('admin/', admin.site.urls),
    path('index/', views.index)
]

api/views.py

from django.http import HttpResponse

# Create your views here.
def index(request):
    return HttpResponse("Hello, world.")

そしたら

$ python manage.py runserver

でサーバーを立ち上げ、http://127.0.0.1:8000/index/にアクセスしHello, world.が表示されれば準備OKです。

それでは非同期メール機能の作成に移ります。

api/mail.pyを新規作成してください。中身は次のように書いてください。

import threading
import queue
import time

q = queue.Queue()

def worker():
    while True:
        item = q.get()
        print(f'Working on {item}')
        time.sleep(3)  # ここで任意の処理を行う
        print(f'Finished {item}')
        q.task_done()

# 処理を並列化するためにスレッドを5つ立てている
for _ in range(5): 
    threading.Thread(target=worker, daemon=True).start()

def add(item):
    q.put(item)

次にviews.pyを以下のように修正します。

from django.http import HttpResponse
from api.mail import add as mail_add

# Create your views here.
def index(request):
    item = request.GET.get("param")
    mail_add(item)
    return HttpResponse("Hello, world.")

この状態でサーバーを再起動し、次のようにリクエストを連続で送信してみます。

curl http://127.0.0.1:8000/index/?param=1
curl http://127.0.0.1:8000/index/?param=2
curl http://127.0.0.1:8000/index/?param=3
curl http://127.0.0.1:8000/index/?param=4
curl http://127.0.0.1:8000/index/?param=5

出力

Working on 1
[15/May/2021 12:12:28] "GET /index/?param=1 HTTP/1.1" 200 13
Working on 2
[15/May/2021 12:12:28] "GET /index/?param=2 HTTP/1.1" 200 13
Working on 3
[15/May/2021 12:12:28] "GET /index/?param=3 HTTP/1.1" 200 13
Working on 4
[15/May/2021 12:12:28] "GET /index/?param=4 HTTP/1.1" 200 13
Working on 5
[15/May/2021 12:12:28] "GET /index/?param=5 HTTP/1.1" 200 13
Finished 1
Finished 2
Finished 3
Finished 4
Finished 5

Finishedがほぼ同時に出力されると思います。time.sleep(3)の部分をメール処理に書き換えれば非同期なメール送信機能の完成です。

【節税】フリーランス一年目ならこれだけやれば80点

こんにちは。フリーランスエンジニア4年目のどすえです。

フリーランス1年目のときに確定申告で苦労したこをと思い出しながら、フリーランス1年目の人が税金・保険料対策として最低限やるべきことをまとめました。

フリーランスはノーガードで確定申告を行うと、これでもかというくらい税金・保険料を支払うことになります。したがって節税が必須です。

節税では大きな部分を間違えなければ致命傷を避けられます。
ここでは手続きをできる限り最低限にして、フリーランス1年目の節税を紹介しようと思います。

ポイントは次の3点です。

  • 大きな支出を経費で精算する。(家賃、水道光熱費、通信費など)
  • 日々の支出はできる限りクレジットカードで支払う。(ただしデビッドカードはNG、理由は後述。)
  • クラウド会計ソフトに加入して確定申告を行う。

さらに余裕があれば

といった流れになります。ここまでできればパイの大きな節税は網羅しているはずです。

目次

確定申告は必ず行う

確定申告とは、1月1日から12月31日までの1年間の所得を合算し、それに対する税額を計算したのち、翌年3月中旬ごろまでに役所に申告するイベントです。

フリーランスなら必ず確定申告を行います。「そろそろ確定申告してくださいね。」と誰かに催促されるわけではないので、うっかりしているとスルーしてしまいそうですが忘れずに行いましょう。

第一ステップとして「確定申告は必ずやる」という認識を持ってください。

税額はどうやって決まるか

税金・保険料は所得に応じて決定されます。所得とは

所得 = 収入 - 控除 - 経費

で計算される金額です。控除は融通が効きづらいですが、経費は自分の支出状況によって変わります。したがって経費を最大限活用することが節税の基本戦略となります。

 

税金・保険料の種類

全てのフリーランスが支払うべき基本の税金・保険料は以下4つです。

まずは、経費を全く計上しなかった場合、これらの税金・保険料がいくらくらい請求されるか手早く計算してみましょう。

所得税

所得税は個人の所得にかかる税金です。所得額が多くなるほど税率が上がる、累進課税制度が適用されます。国税庁HPに税率の試算表が載っています。

f:id:dosuex:20210512011319p:plain
[出典]国税庁HP
事業者によっては業務委託契約であっても源泉徴収を実施しているケースがあります。源泉徴収で徴収された分は所得税として事業者が納めているため、最終的に自身の所得が決定した際は、実際に支払うべき所得税と、源泉徴収で支払い済みの所得税の差額を払う(あるいは還付を受ける)ことになります。

住民税

住民税は自身の住民票のある自治体に払う地方税です。住民税も所得に応じて決定されます。税率はだいたい10%と考えてください。

国民健康保険

フリーランス国民健康保険に加入し、住んでいる自治体に保険料を支払います。保険料は所得、世帯人数、年齢などによって決定されるため、保険料の試算は自治体が作成している保険料の計算表を活用するのが便利です。
保険料は所得に割合をかけて決定される所得割と、一律で固定額が付加される均等割の合計で決まります。 大まかな計算として、所得割が所得の10%、均等割で5万円が付加されると考えてください。

国民年金

フリーランスは厚生年金などに加入できないため、国民年金に加入することになります。
保険料は所得に関わらず一定で、令和3年度分の保険料は月額16,610円です。年間で199,320円です。

試算

私の1年間の収入が300万円だったと仮定して、白色申告で確定申告を行なったとしましょう。

収入=300万
控除=基礎控除38万 + 白色申告控除10万
所得 = 収入 - 控除 = 252万

所得税 住民税 国民健康保険 国民年金
税率 10% 10% 10% + 5万 定額
25万 25万 25万+5万 20万 95万

合計で95万円を支払う試算になります。これでは年の生活費を200万円で慎ましく暮らしたとしても、ギリギリの水準です。

ここで登場するのが経費です。

経費

経費とは事業を行うために支払った費用のことです。経費として計上した支出は、所得額から引くことができるので、所得が減り課税額が減ることに繋がります。

経費の種類

経費として仕分ける際は勘定科目を設定して分類していきます。フリーランスがよく使うであろう勘定科目はこの辺りでしょうか。

勘定科目
旅費交通費 取引先への移動でかかった費用
新聞図書費 技術書などの業務に必要な書籍代
交際費 取引先とのランチ代
消耗品費 ペンやノート
通信費 インターネット回線の利用料金
地代家賃費 オフィスの家賃
水道光熱費 オフィスの水道光熱費

大きい経費は確実に計上する

上記勘定科目の中で割合が大きいのは地代家賃費です。

賃貸を自宅兼オフィスとして使用しているフリーランスの方も多いでしょうから、その場合経費として計上できます。しかし、経費にできるのは事業に関わる部分だけなので、家事按分という考え方が必要になります。

家事按分とは家事、つまり自分が住むために使用している割合と事業で使っている割合で適切に按分するということです。明確に何割と決まっているわけではないので、税務署に按分の根拠を聞かれたとしても合理的な回答ができる割合に設定しましょう。1日のうち事業のために部屋を使っている時間で計算したり、使っている面積で計算したりと様々な按分の考え方があります。家事按分は水道光熱費や通信費にも適用します。

地代家賃費などは計上のハードルが少し高いですが、一気に所得額を下げられるチャンスなので確実に押さえましょう。

例えば、家賃、水道光熱費、通信費などで年間140万円かかっていたとして、家事按分50%で経費計上すると、所得額が70万円減ることになります。この減った金額で元々課税されるはずだった税・保険料がおおざっぱに30%と計算すると、21万円の税金を減らすことができます。

クレカで支払えば領収書は不要!

経費として計上するためには一般的に領主書が必要です。ただし確定申告の際に提出はせず、申告者が手元に保管しておきます。しかし、領収書ベースの経費計上は、領収書本体の物理的な保管、確定申告の際の勘定科目仕分けと、項目ごとの合計金額の計算などを自分で行うことになります。1年間で何百・何千と領収書が発生することを考えるとこれは途方もない作業です。

そこでクレジットカードの出番です。実は領収書の代わりとしてクレジットカードの利用明細が使えます。というのも

  • 店名
  • 日付
  • 商品の内容
  • 金額
  • 購入者の氏名

が記載されていれば領収書の代わりにできるという条件があるためです。これらの条件を満たしていればレシートでもOKです。

カードの利用明細はカード会社に保管されますから、領収書の物理的な保管は不要となります。

ただし、デビットカードで支払った場合は商品の内容が使用明細に記載されないため、領収書代わりに使うことはできません。

クラウド会計でさらに楽に!

さらに、クレジットカードをクラウド会計ソフトと連携させることで、仕分け作業や経費の計算作業も全てクラウド上で完結させることができます。

仕分け作業はスマホのアプリ上で指一本でできます。年度末になると会計ソフトが確定申告書類を全て作ってくれるので、後は印刷して提出すれば終わりです。

この体験は本当に素晴らしく、煩雑で面倒臭い確定申告の作業が極限まで簡略化されています。領収書を紛失したり、計算ミスをしたりという心配もないので、節税のモチベーションも下がりません。まだ利用されていない方は、クラウド会計ソフトを使っていただければ確実に節税額を増やせると思います。おすすめです。

私は過去3年間のクラウド会計ソフト freee を有料会員で利用しています。

仕分けの際に遡れる過去の明細に限りがありますが無料でも使えます。確定申告を行う年初に1ヶ月だけ加入して一気に一年分の取引を登録するという使い方をしている人もいます。いずれにせよ相当な節税効果があります。

また、マネーフォワードクラウドの確定申告サポートも評判が良いです。私は家計簿サービスのほうでマネーフォワードは利用しています。

青色申告

開業届けを役所に提出した個人事業主であれば、青色申告をおこなうことができます。青色申告とは確定申告時に提出する帳簿を複式簿記で記入するなど、少々難易度の高い申告法となっています。それゆえにメリットも大きく、白色申告特別控除が10万円だったのに対し、青色申告特別控除は55万円です。つまり青色申告に切り替えるだけで所得45万円分の追加の節税効果があります。

難易度が高いと言っても、経費の仕分けさえ済んでいれば会計ソフトが自動で作成してくれるので実質難易度は変わりません。開業届を提出済みであれば、青色申告が簡単にできることも、会計ソフト導入のメリットです。

ちなみに私は開業届けもfreeeで作成して提出しました。

まとめ

改めてまとめると、フリーランスが経費を活用して所得を減らし、税・保険料を節税するにはまず以下の3ステップです。

  • 大きな支出を経費で精算する。(家賃、水道光熱費、通信費など)
  • 日々の支出はできる限りクレジットカードで支払う。
  • クラウド会計ソフトに加入して確定申告を行う。

さらに余裕があれば

さいごに

以上、フリーランス一年目ならこれだけやれば80点な節税というテーマでまとめました。この記事を書いたきっかけは、私自身がフリーランス1年目のとき税金・保険料の仕組みが全く分からず、節税がほとんどできぬまま確定申告を行い地獄の2年目を過ごしたからです。

税金関連は法律がダイレクトに関わる話題なので、個別で細かく丁寧に説明してくれる情報源は多いのですが、結局何をしておけば良いのかよくわからなかった思い出があります。

この記事があなたの助けになっていれば嬉しいです。

PythonのpsutilでCPU使用率を並列的に監視する

Pythonで何か処理を走らせている時に, CPUの使用率を定期的に監視するスクリプトを書いた.
マルチスレッドにすることで, メインの処理を止めずに定期実行できる. 重い処理の計算機負荷の時系列データをレポートしたい時とかに使えるかもしれない.
コード全文は最下に.

ライブラリのインストール

psutilというライブラリを使うことでCPUの使用率をPythonで取得することができる. pip経由でインストールできる.

pip install psutil

基本的な使い方は以下の通り.

# percpu=Trueでコアごとの使用率(%)をリストで返す
>>> psutil.cpu_percent(interval=1, percpu=True)
[7.0, 1.0, 8.9, 0.0, 5.0, 1.0, 4.0, 0.0]

# percpu=Falseで平均の使用率となる
>>> psutil.cpu_percent(interval=1, percpu=False)
8.9

# interval(秒)を小さくすると計測時間が短くなるが, 誤差も大きくなる.
>>> psutil.cpu_percent(interval=0.5, percpu=True) 
[6.0, 0.0, 7.8, 1.9, 4.0, 0.0, 5.8, 0.0]

並列処理

今回はネイティブのライブラリであるthreadingを使って並列処理を実装する.
以下のようにスレッドインスタンスの作成と実行を行うことができる.

m = threading.Thread(target=monitor_cpu,args=((initial_time,)))
m.start()

threading.Thread()に渡す引数として, targetには実行するメソッドを渡し, argsには実行メソッドへの引数を渡す.
monitor_cpuが監視用メソッドである.

監視用メソッド

threading.Eventを使って監視用スレッドの管理をする.
threading.Eventの使い方については以下のページが詳しい.
おまいらのthreading.Eventの使い方は間違っている

threading.Eventとはイベントが発生するまでスレッドを待機させ、他のスレッドからイベントを発生させると待機スレッドが再開する、という使い方をする為のクラスです。

上のページでも説明があるように
wait(timeout) : イベントが発生するかtimeout秒経過するまで現在のスレッドを待機させる.
set() : イベントを発生させる.
のふたつを使用する.

def monitor_cpu(initial_time):
    print("START monitor_cpu")
    while not event.wait(1):  # 1秒待機, イベントが発生していなければループ内を実行して再び待機
        elapsed_time = time.time() - initial_time
        cpu_percent = psutil.cpu_percent(percpu=True)
        cpu_percent = '\t'.join(["{:10.4f}".format(v) for v in cpu_percent])
        print("time:", int(elapsed_time), cpu_percent)
    print("END monitor_cpu")  # イベントが発生したらループを抜け実行終了


if __name__=="__main__":
    event = threading.Event()
    initial_time = time.time()
    m = threading.Thread(target=monitor_cpu,args=((initial_time,)))
    m.start()

    # メイン処理

    event.set()  # イベントを発生させる.

monitor_cpu内部では実行間隔を1秒で無限ループを回す. イベントが発生するとループを抜け実行終了する.
ちなみにこの実装はtime.sleep(1)と終了条件を司るフラグを使っても再現できる. こちらのほうがシンプルで良いかもしれない.

def monitor_cpu(initial_time):
    print("START monitor_cpu")
    while flag:
        time.sleep(1)
        elapsed_time = time.time() - initial_time
        cpu_percent = psutil.cpu_percent(percpu=True)
        cpu_percent = '\t'.join(["{:10.4f}".format(v) for v in cpu_percent])
        print("time:", int(elapsed_time), cpu_percent)
    print("END monitor_cpu")

if __name__=="__main__":
    event = threading.Event()
    initial_time = time.time()
    flag = True
    m = threading.Thread(target=monitor_cpu,args=((initial_time,)))
    m.start()
    tmp = 0
    for i in range(100000000):
        tmp = i+i
    flag = False

コード全体

import threading
import time
import psutil

def monitor_cpu(initial_time):
    print("START monitor_cpu")
    while not event.wait(1):
        elapsed_time = time.time() - initial_time
        cpu_percent = psutil.cpu_percent(percpu=True)
        cpu_percent = '\t'.join(["{:10.4f}".format(v) for v in cpu_percent])
        print("time:", int(elapsed_time), cpu_percent)
    print("END monitor_cpu")

if __name__=="__main__":
    event = threading.Event()
    initial_time = time.time()
    m = threading.Thread(target=monitor_cpu,args=((initial_time,)))
    m.start()
    tmp = 0
    for i in range(100000000):
        tmp = i+i
    event.set()

実行結果

左側から時刻(秒), CPU1の使用率(%), CPU2の使用率 ... CPU8の使用率 と並んでいる.

START monitor_cpu
time: 1  52.0     1.0    29.7     1.0    28.7     2.0    30.4     0.0
time: 2  45.0     3.0    36.3     3.0    32.4     1.0    29.4     1.0
time: 3  43.6     2.0    31.0     1.0    33.7     0.0    27.7     2.0
time: 4  45.5     2.0    25.2     1.0    22.8     0.0    19.8     0.0
time: 5  41.6     1.0    26.0     1.0    26.7     1.0    20.8     1.0
time: 6  46.1     5.0    34.7     3.0    38.0     3.0    31.7     4.0
time: 7  63.0    10.0    52.5    10.0    55.4    10.9    55.0    10.0
time: 8  51.5     6.9    36.0     4.9    41.2     5.0    39.6     4.9
time: 9  55.0     0.0    20.6     1.0    26.0     0.0    14.9     0.0
time: 10  49.5    2.9    22.8     1.0    25.2     1.0    23.5     1.0
time: 11  43.6    1.0    32.7     0.0    28.3     0.0    23.8     1.0
time: 12  47.5    2.9    22.8     1.0    20.6     2.0    25.0     2.0
END monitor_cpu