ブログ

割とコンピュータよりの情報をお届けします。

Python

オープンソースのロータダイナミクス解析支援ROSS

ROSSを使うと、参考書片手に手間をかけて計算していた回転体の危険速度の計算などをほぼ実装されたまま利用させていただくことができる。

今まで、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に置き換えていく作業が必要である。(結構面倒)チュートリアルに頼らないのであれば、このバージョン指定インストールは不要である。ただ、一回参考ページに従ってみたい場合には上のようにしないといけない。

≫ Read More

2023/06/21 コンピュータ   TakeMe
Tag:Python

Pythonでsscanfのようなものはないか

「Python sscanf相当」で検索していると「ない」と検索されることが多いけど、parseライブラリにあるような気がする。

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);

≫ Read More

2023/06/21 コンピュータ   TakeMe
Tag:Python

Pythonで3次元プロット と回転 matplotlib 3.7.1

matplotlibの2くらいから「パースペクティブ」の解除は一語に代わっていた。

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()
    

≫ Read More

2023/06/19 コンピュータ   TakeMe
Tag:Python

PyScriptはすぐに使えるのか?

CMSと組み合わせてもPyScriptはすぐに使えるのか非常に疑問だった。
一応バージョン番号付けが1未満になっている意図は察しつつ…。

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

意図通りにはならない。

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

≫ Read More

2023/06/06 コンピュータ   TakeMe
Tag:Python

Python numpyでFFTを実行する.周波数軸はどうしたら4

前の記事ではnumpy.fft.fftfreqをプロットの時の範囲指定に使うことしか考えていなかったが,フィルタなどもできる.

numpy.fft.fftfreqによって生成した周波数軸とnumpy.whereを使うと周波数でフィルタリングが行える.
例えば,
Fzf = np.fft.fft(z);
Fz2f = np.where(np.abs(freq2) > 12, 0, Fzf)
などとすれば,周波数成分を取り除くことができる(0で埋めることができる).
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; # 折り返すのでパワーが2分の1になっている.
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();


≫ Read More

2021/04/03 コンピュータ   TakeMe
Tag:Python

Python numpyでFFTを実行する.周波数軸はどうしたら3

以前の記事「Python numpyでFFTを実行する.周波数軸はどうしたら」「Python numpyでFFTを実行する.周波数軸はどうしたら」ではlinspaceを使って周波数リストを出力していたのをnumpy.fft.fftfreqでもよいと書いていたが,,,

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; # 折り返すのでパワーが2分の1になっている.
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; # 折り返すのでパワーが2分の1になっている.
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();

≫ Read More

2021/03/30 コンピュータ   TakeMe
Tag:Python

Python for .NETの新しい話 Vol. 2

Pythonでsgolayフィルタを使ってみた..NETでもオープンソースのコードがないかなと思っていたのですが,なかなか見つからなかったのでscipyのコードを使うことにしてみた.

自分で探していたけど「SciPy で Savitzky-Golay フィルタ」という参考にmodeなるものも指定できることが分かった.SciPyの方が簡単かな

 

やっぱり簡単.ただし,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;

≫ Read More

2021/02/22 コンピュータ   TakeMe

Python for .NETの新しい話

だいぶん前に「PythonからC#で書いた.NET Framewokのクラスライブラリを読みだす」1つ目2つ目
では.NET Frameworkのプログラムを読み出すことだった。

逆の「C#からPythonの関数を呼び出し」は更新が追い付いていないといって最新版をとってきていた。
多少最新バージョンからバージョンがずれていてもよいならpythonnet_py37_winなどなどnugetで入れられるようになっていた。むしろこちらの方が簡単だな。

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();
                }
            }
        }
    }
}

≫ Read More

2021/02/17 コンピュータ   TakeMe
Tag:Python

WinPython 3.7.4はVSCode付もある

最近はWinPythonの配布ファイルの名前の付け方がわからなくてZeroを選ぶようにしていた.

ところが,最近の配布ファイルを見ているといつのころからなのか不明だが,Visual Studio Codeのポータブル版が同梱されたファイルがある.

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

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

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

≫ Read More

2019/09/29 コンピュータ   TakeMe
Tag:Python

Scipyのinterpolateはデータの外側へも広げられる

Scipyのinterpolateはデータの外側へも広げられることを知った.
fminなどを使用していると,引数が意図しない範囲に行って補間が失敗して止まることがあるのだが,未然に防ぐことができる.

scipyのinterpolateは今まで使ったことがなかったが,事前にx, yの2つの数値列を与えておくと補間するための関数を返してくれるという仕様になっている.この関数として返すという仕様が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次微分係数,...

≫ Read More

2019/07/27 コンピュータ   TakeMe
Tag:Python