ブログ

今まで、GNU OctaveやPythonでプログラムを作って行列計算などをしていたが、ROSSなるものが無料で配布されている。
ただし参考ページ（【ターボ機械シリーズ④】　オープンソースローターダイナミクス解析ソフトROSS）がページ内で参照しているTurboToolsのチュートリアルはross-rotordynamicsのVersion 1.0.0、Numpy 1.2.0以下用になっていてそのまま実行できない。
Numpyの方は軽微な修正でnumpy.complexがなくなったというものなので、complexに置き換えるだけでよい。結構有名なパッケージの仕様変更なので有名である。
ただし、ross-rotordynamicsの方は有名パッケージではないのでversion 1.0.1にバージョン指定をしてインストールした方が良い。GitHubの変更履歴を追ってみると、この後のバージョンで.kxx.plot()の定義がなくなる。

pip install ross-rotordynamics==1.0.1

もしかすると、使用する描画機能にplotlyの出力を抑制するように変化が起こっているかもしれない。
インストールしたrossの中で、numpy.complexを使用している部分を順次complexに置き換えていく作業が必要である。（結構面倒）チュートリアルに頼らないのであれば、このバージョン指定インストールは不要である。ただ、一回参考ページに従ってみたい場合には上のようにしないといけない。

PythonにC言語のsscanf相当のものはないのかということを検索していると、「splitで代用する」ようなものが検索喧嘩の先頭にきていやだ。
標準では厳しいが、parseライブラリにこの辺りの処理の対応がある。
例えば次のような例で使える。

import parse
import numpy as np

if __name__ == '__main__':
    test_string = '2023-06-19T01:20:10.010';
    parsed = parse.parse('{}-{}-{}T{}:{}:{}', '2023-06-19T01:20:10.010');
    
    print(parsed);
    
    parsed = np.array(list(parsed), dtype=float);
    print(parsed);

matplotlibでパースペクティブの解除について、昔は少しコードを書かないといけなかった。
がそれは今は昔の話で、ax.set_proj_type('ortho');で一発となっていた。
まえに記事を書いてから5年近くたっているから当然かな。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def rotX(X):
    theta = X / 180.0 * np.pi;
    Rx = np.array([[1,0,0,0],
                   [0,np.cos(theta),np.sin(theta),0],
                   [0,-np.sin(theta),np.cos(theta),0],
                   [0,0,0,1]]);
    return Rx;

def rotY(Y):
    theta = Y / 180.0 * np.pi;
    Ry = np.array([[np.cos(theta),0,-np.sin(theta),0],
                   [0,1,0,0],
                   [np.sin(theta),0,np.cos(theta),0],
                   [0,0,0,1]]);
    return Ry;

def rotZ(Z):
    theta = Z / 180.0 * np.pi;
    Rz = np.array([[np.cos(theta),np.sin(theta),0,0],
                   [-np.sin(theta),np.cos(theta),0,0],
                   [0,0,1,0],
                   [0,0,0,1]]);
    return Rz;

def trans(x):
    Tr = np.array([[1,0,0,x[0]],
                   [0,1,0,x[1]],
                   [0,0,1,x[2]],
                   [0,0,0,1]]);
    return Tr;

def main():

    r = 50.0;
    theta = np.arange(0.0, 360.0, 30.0);
    x = np.array([]);
    y = np.array([]);
    z = np.array([]);

    for i in range(10):
        x_ = r * np.cos(theta / 180.0 * np.pi);
        y_ = r * np.sin(theta / 180.0 * np.pi);
        z_ = 10.0 * i * np.ones(x_.shape[0]);
        
        x = np.append(x, x_);
        y = np.append(y, y_);
        z = np.append(z, z_);
    
    X = np.array([x, y, z, np.ones(x.shape[0])]);
    X = np.dot(rotX(45.0), X);
    X = np.dot(rotY(45.0), X);
    X = np.dot(rotZ(10.0), X);
    X = np.dot(trans([10, 10, -10]), X);
    
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.set_proj_type('ortho');
    ax.scatter3D(X[0], X[1], X[2]);
    ax.set_aspect('equal')
    plt.grid()
    plt.savefig('plot201803.svg');
    plt.show();

if __name__ == '__main__':
    main()
    

しかも、バックエンドとしてWebAggを使用するとWebブラウザから3D図をグリグリ操作できる(tornadoパッケージも必要)。私はWebView2と組み合わせて.NETアプリとPythonアプリの組み合わせで使ってみようかと思っています。
（バックエンドをWebAggに切り替えただけでは、show()のところでWebブラウザが起動してしまうところが問題だがほかのアプリへの組み込み方があるらしい）

if __name__ == '__main__':
    plt.rcParams['webagg.open_in_browser'] = False;
    plt.rcParams['webagg.address'] = '0.0.0.0';
    plt.switch_backend('WebAgg');
    main()
    

デモのhello worldを張り付けてみたら以下のようになる。最初はpyscriptのライブラリを見込む必要がある。結構ロードするファイルサイズがでかいことが分かったので公開を中止した。
改めてspimle clockの方をインラインフレームで作り直した。

意図通りにはならない。

HTMLではコード中の改行を気にしない仕組みだからでないかと思われる。それ用に作られているCMSのエディタでは改行を気にしない。Pythonは改行とインデントを気にする。（違うのかな）
少し検討が必要かもしれないが、サンプルではMatplotlibの使用例を含んでおり、データをお客に見せるときにこれでよいのであれば使えるかもしれない。
念のため確認だが、依然としてChrome系統のブラウザで動作確認をして活発に更新されているのでわずかな変化で動作に問題が生じる可能性がある。
デモページの3d表示のサンプルも直前まで動いていたような気がしたが、今日は動いていない

numpy.fft.fftfreqによって生成した周波数軸とnumpy.whereを使うと周波数でフィルタリングが行える．
例えば，
Fzf = np.fft.fft(z);
Fz2f = np.where(np.abs(freq2) > 12, 0, Fzf)
などとすれば，周波数成分を取り除くことができる（０で埋めることができる）．
np.fft.ifftで時間領域の信号に戻せば0で埋めた成分は消えている．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 1000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fzf = np.fft.fft(z); # Fzfはfftした結果
Fz = Fzf / z.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);


freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);# 周波数軸の生成

# abs(freq2)が12以上の成分をカット
Fz2f = np.where(np.abs(freq2) > 12, 0, Fzf) # filter
z_filtered = np.fft.ifft(Fz2f);
Fz2 = Fz2f / z.shape[0] * 2;
Fz2[0] = Fz2[0] / 2;
Fz2_abs = np.abs(Fz2);


plt.figure(1)
# 比較FFTプロット
plt.subplot(211);
plotrange = np.where(freq2 >= 0)
plt.plot(freq2[plotrange], Fz_abs[plotrange]);
plt.plot(freq2[plotrange], Fz2_abs[plotrange]);

# 比較波形プロット
plt.subplot(212);
plt.plot(t, z);
plt.plot(t, z_filtered);

plt.savefig("valuestar.work.177.png")

plt.show();

numpy.fft.fftfreqで作る次のようなコードではfreq2は途中でいきなり負に飛ぶことになり，plt.plot(freq2, Fz_abs);とするとプロットが汚くなる．xlimで範囲制限をしても枠外の点とつながった線が気になる．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 1000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp - 1.0/(t_len + t[1] - t[0]), num=len(t));
freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);

plt.figure(1)
plt.plot(freq, Fz_abs);

plt.figure(2)
plt.plot(freq2, Fz_abs);
plt.xlim([0, 125])


plt.show();

たとえば，np.where(freq2 >= 0, freq2, np.nan)で選別してもよい．要素がNaNになっているとプロットに失敗するので描画されない（plt.figure(3)に示した）．
もしくは，np.where(freq2 >=0)でインデックス列の状態にしておいてそれでプロット範囲を周波数で制限することも可能でplt.figure(4)に示した．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 1000);
z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi);

# サンプリング周波数
fsmp = 1 / (t[1] - t[0]);
# 解析時間
t_len = t.max() - t.min();

Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp - 1.0/(t_len + t[1] - t[0]), num=len(t));
freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);

plt.figure(1)
plt.plot(freq, Fz_abs);

plt.figure(2)
plt.plot(freq2, Fz_abs);
plt.xlim([0, 125])

freq3 = np.where(freq2 >= 0, freq2, np.nan)
plt.figure(3)
plt.plot(freq3, Fz_abs);

plt.figure(4)
plotrange = np.where(freq2 >= 0)
plt.plot(freq2[plotrange], Fz_abs[plotrange]);

plt.show();

やっぱり簡単．ただし，Scipyはかなりストレージを使う．これを小さくできたらよいのだが...

C#側のコードは次のような感じにしてみた．

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.IO;
using System.Linq;
using System.Reflection;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using Python.Runtime;

namespace SgolaySample
{
    public partial class Form1 : Form
    {

        public Form1()
        {
            InitializeComponent();

            System.Windows.Forms.DataVisualization.Charting.ChartArea chartArea1 = new System.Windows.Forms.DataVisualization.Charting.ChartArea();
            System.Windows.Forms.DataVisualization.Charting.Legend legend1 = new System.Windows.Forms.DataVisualization.Charting.Legend();
            System.Windows.Forms.DataVisualization.Charting.Series series1 = new System.Windows.Forms.DataVisualization.Charting.Series();
            System.Windows.Forms.DataVisualization.Charting.Series series2 = new System.Windows.Forms.DataVisualization.Charting.Series();
            this.chart1 = new System.Windows.Forms.DataVisualization.Charting.Chart();
            // 
            // chart1
            // 
            chartArea1.Name = "ChartArea1";
            this.chart1.ChartAreas.Add(chartArea1);
            legend1.Name = "Legend1";
            this.chart1.Legends.Add(legend1);
            this.chart1.Location = new System.Drawing.Point(12, 12);
            this.chart1.Name = "chart1";
            series1.ChartArea = "ChartArea1";
            series1.ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Point;
            series1.Legend = "Legend1";
            series1.Name = "Series1";
            series2.ChartArea = "ChartArea1";
            series2.ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Line;
            series2.Legend = "Legend1";
            series2.Name = "Series2";
            this.chart1.Series.Add(series1);
            this.chart1.Series.Add(series2);
            this.chart1.Size = new System.Drawing.Size(446, 206);
            this.chart1.TabIndex = 0;
            this.chart1.Text = "chart1";


            this.Controls.Add(this.chart1);

            string strPath = Environment.GetEnvironmentVariable("PATH");

            string appDir = Directory.GetParent(Assembly.GetExecutingAssembly().Location).FullName + @"\python";

            Environment.SetEnvironmentVariable("PATH", Path.PathSeparator + appDir, EnvironmentVariableTarget.Process);


            List<double> x = new List<double>(); // Python listへの変換予定
            List<double> y = new List<double>(); // Python listへの変換予定

            Random random = new Random(0);
            for (int i = 0; i < 400; i++)
            {
                x.Add(0.01 * i);
                y.Add(Math.Sin(2.0 * Math.PI * x[i] / 1.0) + random.NextDouble());
            }

            using (Py.GIL())
            {
                dynamic sys = Py.Import("sys");
                // sample.pyを置くフォルダをパスに追加
                sys.path.append(Directory.GetParent(Assembly.GetExecutingAssembly().Location).FullName);
            }

            using (Py.GIL()) 
            {
                dynamic sample = Py.Import("sample");
                dynamic ya = sample.filterA(y); // numpy.arrayで返ってくる予定

                for (int i = 0; i < x.Count; i++) {
                    chart1.Series[0].Points.AddXY(x[i], y[i]);
                    chart1.Series[1].Points.AddXY(x[i], ((double[])ya)[i]);
                }
            }
        }
    }
}

sample.pyの例

import numpy as np
from scipy import signal

def filterA(y):
    ya = signal.savgol_filter(y, 101, 5);
    return ya;

nugetパッケージは.NET Frameworkならpythonnet_py37_winだし.NET Coreならpythonnet_netstandard_py37_win
やっぱりnugetパッケージは簡単．ただし，embeddable pythonをビルドターゲットフォルダに持ってきている。この参考でscipyまたはnumpyをインストールした方が便利。

ターゲットに入れて一緒に持ち運び....（scipyまで入れると相当　大きい）

using Python.Runtime;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Reflection;
using System.Text;
using System.Threading.Tasks;

namespace PythonEmbedSample
{
    class Program
    {
        static void Main(string[] args)
        {
            string strPath = Environment.GetEnvironmentVariable("PATH");

            string appDir = Directory.GetParent(Assembly.GetExecutingAssembly().Location).FullName + @"\python";

            Environment.SetEnvironmentVariable("PATH", Path.PathSeparator + appDir, EnvironmentVariableTarget.Process);

            using (Py.GIL())
            {
                PythonEngine.RunSimpleString("print('SAMPLE')");
            }
        }
    }
}

python37._pthの#import siteのコメントアウトを解除して，https://bootstrap.pypa.io/get-pip.pyをもらってきて，...
embeddable pythonのディレクトリでpython get-pip.pyを実行して，
python -m pip install scipyを実行して... scipyを使えるようになった

例えば，sample.pyを作る splineで補間する例を追加してみた．

import numpy as np
from scipy import interpolate

def spline(x, y, xi):
    x = np.array(x);
    y = np.array(y);
    xi = np.array(xi);
    f = interpolate.interp1d(x, y, kind="cubic", fill_value="extrapolate")
    yi = f(xi);
    return yi;

そしてC#の側を次のようにいじってsample.pyのモジュールを読み込むようにする．

using Python.Runtime;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Reflection;
using System.Text;
using System.Threading.Tasks;

namespace PythonEmbedSample
{
    class Program
    {
        static void Main(string[] args)
        {
            string strPath = Environment.GetEnvironmentVariable("PATH");

            string appDir = Directory.GetParent(Assembly.GetExecutingAssembly().Location).FullName + @"\python";

            Environment.SetEnvironmentVariable("PATH", Path.PathSeparator + appDir, EnvironmentVariableTarget.Process);

            using (Py.GIL())
            {
                dynamic sys = Py.Import("sys");
                // sample.pyを置くフォルダをパスに追加
                sys.path.append(Directory.GetParent(Assembly.GetExecutingAssembly().Location).FullName);


                dynamic sample = Py.Import("sample");
                
                List<double> x = new List<double>(); // Python listへの変換予定
                List<double> y = new List<double>(); // Python listへの変換予定

                for (int i = 0; i < 10; i++)
                {
                    x.Add(0.1 * i);
                    y.Add(Math.Sin(2.0 * Math.PI * x[i] / 1.0));
                }

                List<double> xi = new List<double>(); // Python listへの変換予定
                for (int i = 0; i < 100; i++)
                {
                    xi.Add(0.01 * i);
                }
                dynamic yi_ = sample.spline(x, y, xi); // numpy.arrayで返ってくる予定
                List<double> yi = new List<double>(); // List<double>への変換は暗黙的にできないため
                yi.AddRange((double[])yi_);

                {
                    StreamWriter sw = new StreamWriter("sample_init.csv");
                    for (int i = 0; i < x.Count; i++)
                    {
                        sw.WriteLine(string.Format("{0},{1}", x[i], y[i]));
                    }
                    sw.Close();
                }

                {
                    StreamWriter sw = new StreamWriter("sample_interp.csv");
                    for (int i = 0; i < xi.Count; i++)
                    {
                        sw.WriteLine(string.Format("{0},{1}", xi[i], yi[i]));
                    }
                    sw.Close();
                }
            }
        }
    }
}

このVisual Studio CodeはWinPythonと同様USBメモリなどリムーバブルメディアから起動できるようになっている．
Pythonのデバッグ設定関係もセットされていた．（最初の利用時には構成の追加で選択が必要）
Visual Studio Code日本語が必要なら，「Japanese Language Pack for Visual Studio Code」を追加すればよい．

PythonのデバッグにIDLEを使用していたが，ブレイクポイントの設定などでVSCodeの方がよく使用するIDEに似ていて便利かもしれない．

ただし，WinPythonの配布ファイルサイズは結構大きくなってきている．ファイルを展開して準備するまでの時間は相当要する．

scipyのinterpolateは今まで使ったことがなかったが，事前にx, yの２つの数値列を与えておくと補間するための関数を返してくれるという仕様になっている．この関数として返すという仕様がfminなどで使うときの便利さにつながる．
さらに，interpolateはデータの範囲が超えたときにエラーにしない仕様がある．例えば以下のような例である

# -*- coding: utf-8 -*-
import numpy as np
import scipy as scipy
from scipy.optimize import fmin
from scipy import interpolate
import matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0, 100, num=1001, endpoint=True);
y = np.zeros(x.shape);

for i in range(y.shape[0]):
    if i < 500 and i > 50:
        y[i] = 1;
f4 = scipy.interpolate.interp1d(x, y, 'cubic', fill_value='extrapolate')
xi = np.linspace(-100, 200, num=121)

plt.plot(x, y, linewidth='3');
plt.plot(xi, f4(xi))
plt.grid(True)

plt.savefig('20190727test.png', format='png', dpi=100)
plt.show();

さらに数値最小化の追加の例までついかすると例えば以下のようになる．

# -*- coding: utf-8 -*-
import numpy as np
import scipy as scipy
from scipy.optimize import fmin
from scipy import interpolate
import matplotlib
import matplotlib.pyplot as plt

x = np.linspace(0, 100, num=1001, endpoint=True);
y = (x - 40) ** 2

f2 = scipy.interpolate.interp1d(x, y, 'cubic', fill_value='extrapolate')

minimum = scipy.optimize.fmin(f2, 1)

xi = np.linspace(-100, 200, num=121)

plt.plot(x, y, linewidth='3');
plt.plot(xi, f2(xi))
plt.plot(minimum, f2(minimum), '+')
plt.grid(True)

plt.savefig('20190727test2.png', format='png', dpi=100)
plt.show();

導関数を使わない数値最小化ならこれらでもよいが，微分係数を使って勾配を下に下らせるように最適化を目指すアルゴリズムを使用したい場合には，計算が増えてしまう．そこで，UnivariateSplineというものもある．まず上の置き換えだけなら，

from scipy.interpolate import UnivariateSplie
f2 = scipy.interpolate.UnivariateSpline(x, y, s=0, ext=0)

にするだけである．さらにf2.derivatives()なるものが使えるようになる
f2.derivatives(minimum)は[ほぼ0, ほぼ0, 2, ほぼ0]になるはず．順にf2の補間結果，1次微分係数，2次微分係数,...