偶数か奇数か判定するプログラム

言語習得の初期によく出てくる偶数か奇数かを判定するプログラムを改めて。

プログラム

python

pythonevenoroddjudge.py

def is_even(num):
    return num % 2 == 0


num = float(input())

print(is_even(num))

julia

juliaevenoroddjudge.jl

function is_even(num)
    return num % 2 == 0
end
function asknumber()
    print("Enter a number: ")
    parse(Float64, readline())
end
num = asknumber()
is_even(num)

juliaのコードではわざわざfunctionを使っているがif文でも書ける。それについてはまた違う記事で。ほかの言語も後で書きます。

動機

機械学習関連の知識を再確認していて、まずは最も基本的な手法、k-NN法（k最近傍法）を勉強がてら自前で実装してみました。

k-NN法とは

まずはk-NNの概要からおさらいします。とてもシンプルでわかりやすいアルゴリズムなのですぐ理解できると思います。
このアルゴリズムは下の図が全てです。

注目データは中心の星印です。まず与えられたデータに対して全てのデータ間の距離を計算します。例えばk=3のとき、注目データに最も近い上位３個のデータのラベルを見ます。この場合、クラス１が２個、クラス２は１個なので最頻値を取るクラス１に分類されます。一方で、k=6としてみると、クラス１が２個、クラス２が６個となるので、このデータはクラス２に分類されます。
このように与えられたデータは近いデータのラベルの多数決で決定されることになります。
ちなみにk=1のときは最近傍法とも呼ばます。
似たような名前の手法にk-means法がありますが、これは教師なしのクラスタリングの手法で、教師付き分類手法であるk-NN法とは異なった手法です。

実装

それでは実際にPythonで実装していきたいと思います。

データセット

データセットは有名なIris Data Setを使います。
ますiris.dataをダウンロードして実際にデータを読み込みます。

knn.py

data = np.loadtxt('iris.data', delimiter=',', dtype='float64',usecols=[0, 1, 2, 3])
labels = np.loadtxt('iris.data', delimiter=',', dtype='str',usecols=[4])
print(data.shape)
print(labels.shape)

(150, 4)
(150,)

このデータセットは４つの特徴量を持つ１５０個のデータとそのラベル（３クラス）で構成されています。

最近傍探索

今回はこれらのライブラリを使います。

knn.py

import numpy as np
from sklearn import model_selection
from scipy.spatial import distance
from sklearn.metrics import accuracy_score
from statistics import mean, median,variance

まずは与えられた点に対して各データとの距離を計算します。

knn.py

distance_matrix = distance.cdist(test_data, train_data)
indexes = np.argsort(distance_matrix, axis=1)

次に与えられたデータの近傍k個のデータに対して、どのクラスのデータが何個あるのか数えます。

knn.py

class_dict = {}
for label in labels:
    class_dict[label] = 0
class_dict_list = []
for data_num, index in enumerate(indexes[:,:self.k]):
    class_dict_list.append(class_dict.copy())
    for i in index:
        class_dict_list[data_num][self._labels[i]] += 1

最後に最も多かったクラスのラベルを特定します。

knn.py

predict_class = []
for d in class_dict_list:
    max_class = max(d, key=d.get)
    predict_class.append(max_class)

k-NNアルゴリズム自体はこれで全部です。

実行

今回は各クラスをランダムに半分に分けて、学習、実行の流れを２０回繰り返します。実際の実行結果は以下のようになります。（ランダムにデータを分けるので、精度に多少の差は出ます。）

training number 1 ...
knn accuracy : 0.9466666666666667
training number 2 ...
knn accuracy : 0.9333333333333333
training number 3 ...
knn accuracy : 0.9466666666666667
training number 4 ...
knn accuracy : 0.9466666666666667
training number 5 ...
knn accuracy : 0.9333333333333333
training number 6 ...
knn accuracy : 0.92
training number 7 ...
knn accuracy : 0.9466666666666667
training number 8 ...
knn accuracy : 0.9466666666666667
training number 9 ...
knn accuracy : 0.8933333333333333
training number 10 ...
knn accuracy : 0.9466666666666667
training number 11 ...
knn accuracy : 0.96
training number 12 ...
knn accuracy : 0.96
training number 13 ...
knn accuracy : 0.96
training number 14 ...
knn accuracy : 0.96
training number 15 ...
knn accuracy : 0.92
training number 16 ...
knn accuracy : 0.96
training number 17 ...
knn accuracy : 0.92
training number 18 ...
knn accuracy : 0.9866666666666667
training number 19 ...
knn accuracy : 0.9333333333333333
training number 20 ...
knn accuracy : 0.96
=================================================
knn accuracy mean : 0.944
knn accuracy variance : 0.00042292397660818664

このように、高い精度が実現できているのがわかると思います。

感想

シンプルなアルゴリズムだけあって非常に簡単に実装できました。今回のデータのように次元が小さく、少ないデータ数ならば非常に効果的です。
しかし、実際には高次元データに対して次元の呪いの効果でユークリッド距離が意味をなさなくなったり、大規模データセットでは時間計算量、空間計算量ともに実用的ではありません。
このような問題にはベクトルをハッシュ化したり、量子化したりする手法によって、メモリ消費量や計算量を削減する方法が提案されています。
実際のソースコードはgithub上で確認できます。
https://github.com/kotaYkw/machine_learning

4.2 Creating a SeqRecord

4.1へ

Using a SeqRecord object is not very complicated, since all of the information is presented as attributes of the class.
SeqRecord objectの使い方はそんなに複雑ではない、すべての情報はクラスの属性として明示されたからです。

Usually you won’t create a SeqRecord “by hand”, but instead use Bio.SeqIO to read in a sequence file for you (see Chapter 5 and the examples below). However, creating SeqRecord can be quite simple.
通常SeqRecordオブジェクトを自分で生成しない、代わりにBio.SeqIOを使ってシーケンスファイルから読み込みます。（chap5と下のサンプルを参考してください）
しかしながら、SeqRecordの生成はとても簡単です。

4.2.1 SeqRecord objects from scratch

To create a SeqRecord at a minimum you just need a Seq object:
Seqクラスを使うだけで最小限のSeqRecordオブジェクトを作れます。

>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq)

Additionally, you can also pass the id, name and description to the initialization function, but if not they will be set as strings indicating they are unknown, and can be modified subsequently:
さらに、初期化の時にid、名前と解説を代入することができる。もしそうしなければ、それらの値はunknownで初期化されるが、後で修正しても大丈夫です。

>>> simple_seq_r.id
'<unknown id>'
>>> simple_seq_r.id = "AC12345"
>>> simple_seq_r.description = "Made up sequence I wish I could write a paper about"
>>> print(simple_seq_r.description)
Made up sequence I wish I could write a paper about
>>> simple_seq_r.seq
Seq('GATC')

Including an identifier is very important if you want to output your SeqRecord to a file.
You would normally include this when creating the object:
もしSeqRecordをファイルとして出力したいなら、識別子を付けることがとても大事になります。通常はオブジェクトを生成する時に付け加えます。

>>> from Bio.Seq import Seq
>>> simple_seq = Seq("GATC")
>>> from Bio.SeqRecord import SeqRecord
>>> simple_seq_r = SeqRecord(simple_seq, id="AC12345")

As mentioned above, the SeqRecord has an dictionary attribute annotations.
This is used for any miscellaneous annotations that doesn’t fit under one of the other more specific attributes.
上で述べたようにSeqRecordはアノテーションという辞書型の属性を持っています。
これは他の明確な属性にはまらない雑多なアノテーションに対して使います。

Adding annotations is easy, and just involves dealing directly with the annotation dictionary:
アノテーションの追加は簡単で、直接辞書の中に追加すればいいです。

>>> simple_seq_r.annotations["evidence"] = "None. I just made it up."
>>> print(simple_seq_r.annotations)
{'evidence': 'None. I just made it up.'}
>>> print(simple_seq_r.annotations["evidence"])
None. I just made it up.

Working with per-letter-annotations is similar, letter_annotations is a dictionary like attribute which will let you assign any Python sequence (i.e. a string, list or tuple) which has the same length as the sequence:
per-letter-annotationsの使い方も似ていて、letter_annotationsは辞書型でシーケンスデータを割り当てることができます。(i.e. 文字列, リスト あるいはタプル)

>>> simple_seq_r.letter_annotations["phred_quality"] = [40, 40, 38, 30]
>>> print(simple_seq_r.letter_annotations)
{'phred_quality': [40, 40, 38, 30]}
>>> print(simple_seq_r.letter_annotations["phred_quality"])
[40, 40, 38, 30]

The dbxrefs and features attributes are just Python lists, and should be used to store strings and SeqFeature objects (discussed later in this chapter) respectively.
データベースへのクロスリファレンス（.dbxrefs）とフィーチャー情報はリスト型、それぞれ文字列およびSeqFeatureオブジェクトを格納するためいに使います。 (この章の後で述べます)

4.2.2 SeqRecord objects from FASTA files

This example uses a fairly large FASTA file containing the whole sequence for Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, originally downloaded from the NCBI.
このサンプルではYersinia pestis菌の全シーケンスが含まれるかなりでかいFASTAファイルを使います。NCBIからダウンロードしました。

This file is included with the Biopython unit tests under the GenBank folder, or online NC_005816.fna from our website.
Biopythonユニットテストを含むファイルはGenBankフォルダにあります、あるいは私たちのサイトからもファイルを取得できます。

The file starts like this - and you can check there is only one record present (i.e. only one line starting with a greater than symbol):
このファイルは下記のような形で始まります - 今は一行しか載せてないですが（i.e.大なり記号から始まる一行）

>gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus ... pPCP1, complete sequence
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...

Back in Chapter 2 you will have seen the function Bio.SeqIO.parse(...) used to loop over all the records in a file as SeqRecord objects.
Chapter 2で見たように関数Bio.SeqIO.parse(...)はSeqRecordオブジェクトファイルの全レコードをループします。

The Bio.SeqIO module has a sister function for use on files which contain just one record which we’ll use here (see Chapter 5 for details):
ここで使うのは1レコードしか持たないファイルを読み込むためのBio.SeqIOモジュールにある類似の姉妹関数です。（詳細は第5章で述べます。）
※意味不明です。

Now, let’s have a look at the key attributes of this SeqRecord individually – starting with the seq attribute which gives you a Seq object:
では、SeqRecordのメイン属性を見てみましょう - seq属性にアクセスしたらSeqオブジェクトを返してくれます。

>>> record.seq
Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG', SingleLetterAlphabet())

Here Bio.SeqIO has defaulted to a generic alphabet, rather than guessing that this is DNA.
Bio.SeqIOが遺伝子アルファベットを返してくれた、推測するまでもなくそれはDNAです。

If you know in advance what kind of sequence your FASTA file contains, you can tell Bio.SeqIO which alphabet to use (see Chapter 5).
もしあなたが予めFASTAファイルにどんなシーケンスデータが含まれると分かったら、Bio.SeqIOモジュールにどのアルファベットを使うのかを指定できます(5章を参照)。

Next, the identifiers and description:
次は「識別子」と「説明」：

>>> record.id
'gi|45478711|ref|NC_005816.1|'
>>> record.name
'gi|45478711|ref|NC_005816.1|'
>>> record.description
'gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus ... pPCP1, complete sequence'

As you can see above, the first word of the FASTA record’s title line (after removing the greater than symbol) is used for both the id and name attributes.
上で見たように、FASTAファイル内のレコード行のタイトルの初文字（大なり記号>を取り除いた後）はidとname属性に両方に使われています。

The whole title line (after removing the greater than symbol) is used for the record description. This is deliberate, partly for backwards compatibility reasons, but it also makes sense if you have a FASTA file like this:
タイトル全行（大なり記号>を取り除いた後）はレコードの説明に使われます。これはバックエンドとの互換性を取るための理由などがあります。しかも以下のようなFASTAファイルの場合はなおさらです。

>Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1
TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGGGGGTAATCTGCTCTCC
...

Note that none of the other annotation attributes get populated when reading a FASTA file:
注意するのはFASTAファイルを投入する初期では、他の属性は空となります。

In this case our example FASTA file was from the NCBI, and they have a fairly well defined set of conventions for formatting their FASTA lines.
このコースでは私たちのサンプルFASTAファイルはNCBIから取得していて、FASTAレコードのフォーマットはかなり整備されています。

This means it would be possible to parse this information and extract the GI number and accession for example.
However, FASTA files from other sources vary, so this isn’t possible in general.
つまり情報の分析およびGIナンバー、付加情報の抽出は可能になります。
しかし他のソース源では難しいことになるでしょう。

4.2.3 SeqRecord objects from GenBank files

As in the previous example, we’re going to look at the whole sequence for Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, originally downloaded from the NCBI, but this time as a GenBank file.
前の例で、NCBIからダウンロードしたYersinia pestis菌の全シーケンスを見ました。今回はGenBankバージョンを見ましょう。

Again, this file is included with the Biopython unit tests under the GenBank folder, or online NC_005816.gbfrom our website.
同じく、GenBankフォルダからBiopythonユニットテスト付きのファイルかを使うか、サイトから直接ダウンロードします。

This file contains a single record (i.e. only one LOCUS line) and starts:
このファイルはシングルレコード(i.e. LOCUS 行はただ１行だけ)およぶstartsを含みます
※starts意味不明

LOCUS       NC_005816               9609 bp    DNA     circular BCT 21-JUL-2008
DEFINITION  Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete
            sequence.
ACCESSION   NC_005816
VERSION     NC_005816.1  GI:45478711
PROJECT     GenomeProject:10638
...

Again, we’ll use Bio.SeqIO to read this file in, and the code is almost identical to that for used above for the FASTA file (see Chapter 5 for details):
同じく、Bio.SeqIOを使ってファイルを読み込み、FASTAファイルの時とほぼ同じやり方です。(詳細は第5章を参照):

>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.gb", "genbank")
>>> record
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG',
IUPACAmbiguousDNA()), id='NC_005816.1', name='NC_005816',
description='Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence.',
dbxrefs=['Project:10638'])

You should be able to spot some differences already! But taking the attributes individually, the sequence string is the same as before, but this time Bio.SeqIO has been able to automatically assign a more specific alphabet (see Chapter 5 for details):
すでにいくつかの違いを発見したはずですが、1個1個見ていきましょう、シーケンスに文字列は前と同じだが、今回Bio.SeqIOはもっと正確にアルファベットの種類を判定できました。（詳細は第5章を参照）

The name comes from the LOCUS line, while the id includes the version suffix. The description comes from the DEFINITION line:
名前はLOCUS行、idがバージョン情報を含んでいます。説明はDEFINITION行から取りました。

GenBank files don’t have any per-letter annotations:
GenBankファイルにper-letterアノテーションは存在しないです。

>>> record.letter_annotations
{}

Most of the annotations information gets recorded in the annotations dictionary, for example:
多くのアノテーション情報はアノテーション辞書に収録されています。たとえば：

>>> len(record.annotations)
11
>>> record.annotations["source"]
'Yersinia pestis biovar Microtus str. 91001'

The dbxrefs list gets populated from any PROJECT or DBLINK lines:
リストdbxrefsはPROJECTあるいはDBLINK行の内容を反映します。

>>> record.dbxrefs
['Project:10638']

Finally, and perhaps most interestingly, all the entries in the features table (e.g. the genes or CDS features) get recorded as SeqFeature objects in the features list.
最後に、一番興味深いのは、フィーチャーテーブルのすべてのエントリ(e.g. genesあるいはCDSフィーチャー)がSeqFeatureオブジェクトとしてフィーチャーリストに保存されます。

何度shutdownしてもゾンビのように生き返る怖いスクリプト

IOS-XEでPythonスクリプトを動作する基本設定はCisco覚書_Pythonでconfig投入参照

スクリプト配置パス

パス：/flash/
ファイル：eem_script.py
guestshellに入ってファイル作成、またはファイルを転送

isp-sw#guestshell
[guestshell@guestshell ~]$ 
[guestshell@guestshell ~]$ cd /flash/
[guestshell@guestshell flash]$ vi eem_script.py 

スクリプトの動作

Loopback55のIFをno shutしてsh ip int briefでIF状態を表示する、という内容
今回はスクリプト実行のトリガーはEEMアプレットで指定する

eem_script.py

import sys
from cli import cli,clip,configure,configurep, execute, executep
intf= sys.argv[1:]
intf = ''.join(intf[0])

print ('This script is going to unshut interface %s and then print show ip interface brief'%intf)

if intf == 'loopback55':
        configurep(["interface loopback55","no shutdown","end"])
else :
        cmd='int %s,no shut ,end' % intf
configurep(cmd.split(','))
executep('show ip interface brief')

手動でスクリプトを起動する方法

EEMアプレットからだけではなく、EXECモードで手動でスクリプトを実行することもできる
※スクリプト中のsh ip int briefコマンドの実行はEXECモードからなのでEEMアプレットでは実行されない。

isp-sw#guestshell run python /flash/eem_script.py loop55
This script is going to unshut interface loop55 and then print show ip interface brief
Line 1 SUCCESS: int loop55
Line 2 SUCCESS: no shut 
Line 3 SUCCESS: end
Interface              IP-Address      OK? Method Status                Protocol
Vlan1                  unassigned      YES NVRAM  administratively down down    
---省略---
Loopback13             10.13.13.13     YES manual up                    up      
Loopback55             10.55.55.55     YES TFTP   up                    up      

EEMアプレット作成

設定モードにて
admindownのsyslogをトリガーにadminUPするスクリプトを実行する内容

isp-sw(config)#event manager applet intshut
isp-sw(config-applet)# event syslog pattern "Interface Loopback55, changed state to administratively down"
isp-sw(config-applet)# action 0.0 cli command "en"
isp-sw(config-applet)# action 1.0 cli command "guestshell run python /flash/eem_script.py loop55"
isp-sw(config-applet)#end

スクリプトのパスをconfig上で指定

設定モードにて
配置パスを指定

isp-sw(config)#event manager directory user policy "flash:/"

EEM起動状態確認

EXECモードにて
正常に登録されていれば表示される

isp-sw#sh event manager policy registered 
No.  Class     Type    Event Type          Trap  Time Registered           Name
1    applet    user    syslog              Off   Thu Mar 28 16:37:52 2020  intshut
 pattern {Interface Loopback55, changed state to administratively down}
 maxrun 20.000
 action 0.0 cli command "en"
 action 1.0 cli command "guestshell run python /flash/eem_script.py loop55"

手動shutdownしてみる

設定モードでshutdownコマンドを投入、ter monでログをターミナルに出力させて動作確認
ログからshurdown後に自動でUPしている。

isp-sw(config-if)#shutdown 
isp-sw(config-if)#
Mar 28 16:56:22.428 JST: %LINEPROTO-5-UPDOWN: Line protocol on Interface Loopback55, changed state to down
Mar 28 16:56:22.428 JST: %LINK-5-CHANGED: Interface Loopback55, changed state to administratively down
isp-sw(config-if)#
Mar 28 16:56:29.880 JST: %LINEPROTO-5-UPDOWN: Line protocol on Interface Loopback55, changed state to up
isp-sw(config-if)#
Mar 28 16:56:29.880 JST: %LINK-3-UPDOWN: Interface Loopback55, changed state to up
isp-sw(config-if)#

イベントが発生した履歴を確認

EXECモードにて
successでEEMにより正常実行している

isp-sw#sh event manager history events 
No.  Job Id Proc Status   Time of Event            Event Type        Name
1    1      Actv success  Thu Mar28 16:36:19 2020  syslog            applet: intshut 
2    2      Actv success  Thu Mar28 16:43:43 2020  syslog            applet: intshut 
3    3      Actv success  Thu Mar28 16:56:22 2020  syslog            applet: intshut 

---

確認機種とバージョン

機種:Catalyst3850
Ver:Cisco IOS XE　16.06.05　Everest

---

参考リンク

Programmability Configuration Guide, Cisco IOS XE Everest 16.6.x

この記事は言語処理100本ノック 2020【第7章: 単語ベクトル】の続きです。

この記事では、第8章(70-79)の機械学習について扱います。

リンク

この記事にはコードのみを載せました。問題文や解き方の補足は下記のリンクを参考にしてください。

言語処理100本ノック2020 第8章: ニューラルネット

第８章: ニューラルネット

70. 単語ベクトルの和による特徴量

import pandas as pd
import gensim
import numpy as np
train = pd.read_csv('train.txt',sep='\t',header=None)
valid = pd.read_csv('valid.txt',sep='\t',header=None)
test = pd.read_csv('test.txt',sep='\t',header=None)
model = gensim.models.KeyedVectors.load_word2vec_format('GoogleNews-vectors-negative300.bin', binary=True)

d = {'b':0, 't':1, 'e':2, 'm':3}
y_train = train.iloc[:,0].replace(d)
y_train.to_csv('y_train.txt',header=False, index=False)
y_valid = valid.iloc[:,0].replace(d)
y_valid.to_csv('y_valid.txt',header=False, index=False)
y_test = test.iloc[:,0].replace(d)
y_test.to_csv('y_test.txt',header=False, index=False)

def write_X(file_name, df):
    with open(file_name,'w') as f:
        for text in df.iloc[:,1]:
            vectors = []
            for word in text.split():
                if word in model.vocab:
                    vectors.append(model[word])
            if (len(vectors)==0):
                vector = np.zeros(300)
            else:
                vectors = np.array(vectors)
                vector = vectors.mean(axis=0)
            vector = vector.astype(np.str).tolist()
            output = ' '.join(vector)+'\n'
            f.write(output)
write_X('X_train.txt', train)
write_X('X_valid.txt', valid)
write_X('X_test.txt', test)

71. 単層ニューラルネットワークによる予測

import torch
import numpy as np
X_train = np.loadtxt(base+'X_train.txt', delimiter=' ')
X_train = torch.tensor(X_train, dtype=torch.float32)
W = torch.randn(300, 4)
softmax = torch.nn.Softmax(dim=1)
print (softmax(torch.matmul(X_train[:1], W)))
print (softmax(torch.matmul(X_train[:4], W)))

72. 損失と勾配の計算

y_train = np.loadtxt(base+'y_train.txt')
y_train = torch.tensor(y_train, dtype=torch.int64)
loss = torch.nn.CrossEntropyLoss()
print (loss(torch.matmul(X_train[:1], W),y_train[:1]))
print (loss(torch.matmul(X_train[:4], W),y_train[:4]))

ans = [] # 以下、確認
for s,i in zip(softmax(torch.matmul(X_train[:4], W)),y_train[:4]):
  ans.append(-np.log(s[i]))
print (np.mean(ans))

73. 確率的勾配降下法による学習

from torch.utils.data import TensorDataset, DataLoader
class LogisticRegression(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 4),
        )
    def forward(self, X):
        return self.net(X)

model = LogisticRegression()
ds = TensorDataset(X_train, y_train)
# DataLoaderを作成
loader = DataLoader(ds, batch_size=1, shuffle=True)

loss_fn = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)

for epoch in range(10):
    for xx, yy in loader:
        y_pred = model(xx)
        loss = loss_fn(y_pred, yy)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

74. 正解率の計測

def accuracy(pred, label):
  pred = np.argmax(pred.data.numpy(), axis=1)
  label = label.data.numpy()
  return (pred == label).mean()


X_valid = np.loadtxt(base+'X_valid.txt', delimiter=' ')
X_valid = torch.tensor(X_valid, dtype=torch.float32)
y_valid = np.loadtxt(base+'y_valid.txt')
y_valid = torch.tensor(y_valid, dtype=torch.int64)

pred = model(X_train)
print (accuracy(pred, y_train))
pred = model(X_valid)
print (accuracy(pred, y_valid))

75. 損失と正解率のプロット

%load_ext tensorboard
!rm -rf ./runs
%tensorboard --logdir ./runs
from torch.utils.tensorboard import SummaryWriter
writer = SummaryWriter()
from torch.utils.data import TensorDataset, DataLoader
class LogisticRegression(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 4),
        )
    def forward(self, X):
        return self.net(X)

model = LogisticRegression()
ds = TensorDataset(X_train, y_train)
# DataLoaderを作成
loader = DataLoader(ds, batch_size=1, shuffle=True)

loss_fn = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)

for epoch in range(10):
    for xx, yy in loader:
        y_pred = model(xx)
        loss = loss_fn(y_pred, yy)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    with torch.no_grad():
      y_pred = model(X_train)
      loss = loss_fn(y_pred, y_train) 
      writer.add_scalar('Loss/train', loss, epoch)
      writer.add_scalar('Accuracy/train', accuracy(y_pred,y_train), epoch)

      y_pred = model(X_valid)
      loss = loss_fn(y_pred, y_valid)
      writer.add_scalar('Loss/valid', loss, epoch)
      writer.add_scalar('Accuracy/valid', accuracy(y_pred,y_valid), epoch)

76. チェックポイント

from torch.utils.data import TensorDataset, DataLoader
class LogisticRegression(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 4),
        )
    def forward(self, X):
        return self.net(X)

model = LogisticRegression()
ds = TensorDataset(X_train, y_train)
# DataLoaderを作成
loader = DataLoader(ds, batch_size=1, shuffle=True)

loss_fn = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)

for epoch in range(10):
    for xx, yy in loader:
        y_pred = model(xx)
        loss = loss_fn(y_pred, yy)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    with torch.no_grad():
      y_pred = model(X_train)
      loss = loss_fn(y_pred, y_train) 
      writer.add_scalar('Loss/train', loss, epoch)
      writer.add_scalar('Accuracy/train', accuracy(y_pred,y_train), epoch)

      y_pred = model(X_valid)
      loss = loss_fn(y_pred, y_valid)
      writer.add_scalar('Loss/valid', loss, epoch)
      writer.add_scalar('Accuracy/valid', accuracy(y_pred,y_valid), epoch)

      torch.save(model.state_dict(), base+'output/'+str(epoch)+'.model')
      torch.save(optimizer.state_dict(), base+'output/'+str(epoch)+'.param')

77. ミニバッチ化

import time
from torch.utils.data import TensorDataset, DataLoader
class LogisticRegression(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 4),
        )
    def forward(self, X):
        return self.net(X)

model = LogisticRegression()
ds = TensorDataset(X_train, y_train)
loss_fn = torch.nn.CrossEntropyLoss()


ls_bs = [2**i for i in range(15)]
ls_time = []
for bs in ls_bs:
  loader = DataLoader(ds, batch_size=bs, shuffle=True)
  optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)
  for epoch in range(1):
      start = time.time()
      for xx, yy in loader:
          y_pred = model(xx)
          loss = loss_fn(y_pred, yy)
          optimizer.zero_grad()
          loss.backward()
          optimizer.step()
  ls_time.append(time.time()-start)
print (ls_time)

78. GPU上での学習

import time
from torch.utils.data import TensorDataset, DataLoader
class LogisticRegression(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 4),
        )
    def forward(self, X):
        return self.net(X)

model = LogisticRegression()
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = model.to(device)

ds = TensorDataset(X_train.to(device), y_train.to(device))
loss_fn = torch.nn.CrossEntropyLoss()

ls_bs = [2**i for i in range(15)]
ls_time = []
for bs in ls_bs:
  loader = DataLoader(ds, batch_size=bs, shuffle=True)
  optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)
  for epoch in range(1):
      start = time.time()
      for xx, yy in loader:
          y_pred = model(xx)
          loss = loss_fn(y_pred, yy)
          optimizer.zero_grad()
          loss.backward()
          optimizer.step()
  ls_time.append(time.time()-start)
print (ls_time)

79. 多層ニューラルネットワーク

import time
from torch.utils.data import TensorDataset, DataLoader
class MLP(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(300, 32),
            torch.nn.ReLU(),
            torch.nn.Linear(32, 4),
        )
    def forward(self, X):
        return self.net(X)

model = MLP()
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = model.to(device)

ds = TensorDataset(X_train.to(device), y_train.to(device))
loss_fn = torch.nn.CrossEntropyLoss()

loader = DataLoader(ds, batch_size=1024, shuffle=True)
optimizer = torch.optim.SGD(model.net.parameters(), lr=1e-1)
for epoch in range(100):
    start = time.time()
    for xx, yy in loader:
        y_pred = model(xx)
        loss = loss_fn(y_pred, yy)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    with torch.no_grad():
      y_pred = model(X_train.to(device))
      loss = loss_fn(y_pred, y_train.to(device)) 
      writer.add_scalar('Loss/train', loss, epoch)
      train_acc = accuracy(y_pred.cpu(),y_train.cpu())
      writer.add_scalar('Accuracy/train', acc, epoch)

      y_pred = model(X_valid.to(device))
      loss = loss_fn(y_pred, y_valid.to(device))
      writer.add_scalar('Loss/valid', loss, epoch)
      valid_acc = accuracy(y_pred.cpu(),y_valid.cpu())
      writer.add_scalar('Accuracy/valid', acc, epoch)
print (train_acc, valid_acc)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import cv2
import time

def conv():
    # 閾値の設定
    threshold = 100

    # 二値化(閾値100を超えた画素を255にする。)
    ret, img_thresh = cv2.threshold(img, threshold, 255, cv2.THRESH_BINARY)

    # 二値化画像の表示
    cv2.imshow("img_th", img_thresh)
    cv2.waitKey()
    cv2.destroyAllWindows()

def fps(video):
    fps = video.get(cv2.CAP_PROP_FPS)
    print("FPSの設定値、video.get(cv2.CAP_PROP_FPS) : {0}".format(fps))

    # 取得するフレームの数
    num_frames = 120

    print("取得中 {0} frames".format(num_frames))

    # 開始時間
    start = time.time()

    # フレームを取得する
    #for i in range(0, num_frames):
    #   ret, frame = video.read()

    # 終了時間
    end = time.time()

    # Time elapsed
    seconds = end - start
    print("経過時間: {0} seconds".format(seconds))

    # Calculate frames per second
    fps = num_frames / seconds
    print("計算したFPS : {0}".format(fps))

class FPS:
    def __init__(self):
        self.flag = False
        self.start = 0
        self.end = 0
        self.fps = 0
        self.framecnt = 0

    def calc_fps(self):
        if self.flag == False:
            self.start = time.time()
            self.flag = True        
        else:
            diff = time.time() - self.start
            if diff > 1.0:
                self.fps = self.framecnt / diff
                self.framecnt = 0
                self.start =  time.time()

        self.framecnt += 1
        return self.fps

def test1(mode):
    URL = "http://172.23.64.38:8081/?action=stream"
    winname ="winname"
    cv2.namedWindow(winname, cv2.WINDOW_NORMAL)

    try :
        if mode == 0:
            s_video = cv2.VideoCapture(0)
        else:
            s_video = cv2.VideoCapture(URL)


        # 閾値の設定
        threshold = 200
        fpsobj = FPS()

        while True:

            start = time.time()
            end = time.time()

            seconds = end - start

            ret, srcimg = s_video.read()
            #cv2.imshow(winname,img)
            key = cv2.waitKey(1) & 0xff
            if key == ord('q'): break

            # 二値化(閾値100を超えた画素を255にする。)
            #cv2.threshold(画像, 閾値, 閾値を超えた場合に変更する値, 二値化の方法)
            ret, img = cv2.threshold(srcimg, threshold, 255, cv2.THRESH_BINARY)

            ksize=7 # アパーチャーサイズ 3, 5, or 7 など 1 より大きい奇数。数値が大きいほどぼかしが出る。
            img = cv2.medianBlur(img,ksize)

            img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # RGB2〜 でなく BGR2〜 を指定
            ret, img = cv2.threshold(img, 20, 255, cv2.THRESH_BINARY)


            #print("{}".format( fpsobj.calc_fps()))
            if ret:
                # 二値化画像の表示
                cv2.imshow(winname, img)
    except:
        pass

    cv2.destroyAllWindows()




if __name__ == "__main__":
    test1(1)

C++

#include <opencv2/core.hpp>
#include <opencv2/videoio.hpp>
#include <opencv2/highgui.hpp>
#include <iostream>
#include <stdio.h>

#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/photo.hpp"

//    http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_thresholding/py_thresholding.html

        Mat dst;
        int criteria = 100;
        int max = 255;


        for ( ;;) {

            cap.read( frame );
            cv::threshold( frame, dst, criteria, max, THRESH_BINARY );       // binalized


            int ksize = 7;// # アパーチャーサイズ 3, 5, or 7 など 1 より大きい奇数。数値が大きいほどぼかしが出る。
            cv::medianBlur( dst, dst, ksize );

            cv::cvtColor( dst, dst, COLOR_BGR2GRAY );// # RGB2〜 でなく BGR2〜 を指定
            cv::threshold( dst, dst, 20, max, THRESH_BINARY );


            imshow( "Live", dst );
            if ( waitKey( 5 ) >= 0 ) {
                break;
            }
        }

もし、あなたの目の前に亡骸があったら？

(引用;名探偵コナン「甘く冷たい宅配便」)
コナン君が遺体を見つけた時、かなり最初の段階で死後経過時間について考えますよね？死亡推定時刻が分かれば大きく容疑者を絞り込むことができますしね＼(^o^)／
そこで今回は死後経過時間について詳しい求め方を紹介します！

例題

室温16℃に設定された冷蔵車に男性の遺体があった。男性の体重は86kg、コナン君が遺体の直腸内温度を計ったところ、27℃であった。男性の死後経過時間を求めよ。

ニュートンの冷却の法則

まず、簡略化した事象を考えてみましょう。コップに入れられた熱いコーヒーが周囲の温度（空気）によって冷まされる問題を考えます。この時、次の「ニュートンの冷却の法則」が近似的に成り立つと言われています。
ニュートンの冷却の法則を表す微分方程式は

$\frac{dT}{dt} = -γ(T - T_s)$ (1)

で表されます。ここで、各文字は
$T$：湯の温度
$T_s$：周囲の温度
$\gamma$：冷却定数
を表しています。この微分法方程式は解析的に解くことができます。まず、(1)式の両辺を$T-T_s$で割ります。
$\frac{1}{T-T_s}dT = -γdt$

両辺を積分して積分定数をCとすると、
$log (T - T_s) = -γt + C$
となります。この式をTについて解くと下式のようになります。
$T = T_s + e^{-γt + C}$ (2)

初期条件$T(0)=T_0$を課すと、積分定数Cは
$C = log (T_0-T_s)$

となります。これを(2)式に代入すれば、解析解は以下で与えられます。
$T(t) = T_s + (T_0 - T_s) e^{-γt}$

（また、これを変形して）
$\frac{T(t)-T_s}{T_0-T_s} = e^{-γt}$

直腸温の連続測定による死後経過時間の推算

次に、ニュートンの冷却の法則を利用して死後経過時間を求めます。「直腸温の連続測定による死後経過時間の推算（早大院先進理工）」によると遺体の場合、微分方程式は次の様に変形されます。
(1) $\frac{T_r-T_e}{T_o-T_e}=1.25e^{B*t}-0.25e^{5B*t}$ $(T_e\le23.2)$
(2) $\frac{T_r-T_e}{T_o-T_e}=1.11e^{B*t}-0.11e^{10B*t}$ $(T_e\ge23.2)$
(3) $B=-1.2815(bW)^{-0.625}+0.0284$
ただし、
$T_r$は直腸温℃
$T_e$は環境温℃
$T_o$は死亡時の直腸温であり37.2℃とする
$T_e$は環境温℃で測定時間の平均値で一定とする
$W$死体の体重kg
$b=1$は補正係数
$B$はNewtonの冷却定数
$t$は死後経過時間(h)
とします。
これを元に早速今回の例題を解きましょう！

import sympy

def death_time(Tr,Te,W):
    Tr=Tr
    Te=Te
    W=W
    B = sympy.Symbol('B')
    t = sympy.Symbol('t')
    To = 37.2
    b = 1
    E = sympy.S.Exp1
    equation1 = 1.25*(E**(B*t))-0.25*(E**(5*B*t))-(Tr-Te)/(To-Te)
    equation2 = 1.11*(E**(B*t))-0.11*(E**(10*B*t))-(Tr-Te)/(To-Te)
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    if Te <= 23.2:
        print(sympy.solve([equation1,equation3]))
    else:
        print(sympy.solve([equation2,equation3]))    
death_time(27,16,86) #今回の問題の設定       

$ $参考文献によると$T_e$は環境温の平均を使いますが、今回は発見時の室温の16℃を$T_e$としました。
結果は言わなくても分かると思いますが、scipyでこんな累乗変数を含む計算を直接できるわけもなく、朝起きても計算を続けていました(泣)。参考文献通りに試行錯誤法の考えの元で解くことにしました。

(引用Krylov Subspace Accelerated Newton Algorithm: Application to Dynamic Progressive Collapse Simulation of Frames)

関数の微分を求め方程式を反復の都度解く処理が大きい為、上図の様に微分を差分で近似させたクリロフ部分交換法を用いて計算することにしました。

import math
from scipy.optimize import newton_krylov
import sympy

#(3)式に相当する。体重(kg)からNewtonの冷却定数Bを算出。To(死亡時直腸温)は37.2とする
def predict_B(W):
    W=W
    B = sympy.Symbol('B')
    b = 1
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    B = sympy.solve([equation3])[B]
    return B  

#(1)(2)式に相当。環境温Teで場合分け
def predict_deathtime(W,Tr,Te,ta):
    Bx = predict_B(W)
    Tr,Te,ta = Tr,Te,ta
    E = sympy.S.Exp1
    To = 37.2
    def F(x):
        if Te <= 23.2:
            return 1.25*(E**(Bx*x))-0.25*(E**(5*Bx*x))-(Tr-Te)/(To-Te)              
        else:
            return 1.11*(E**(Bx*x))-0.11*(E**(10*Bx*x))-(Tr-Te)/(To-Te)
    guess = ta
    sol = newton_krylov(F, guess, method='lgmres')
    print(sol)

#(体重,直腸温,環境温,死後経過時間の初期値)の順で入力
#死後経過時間の初期値を0とするとヤコビアンが求められないので注意する
predict_deathtime(86,27,16,2)

出力結果17.19191007875492

よって、今回の例題では17時間ほど前に死亡したことが分かりました！

参考文献

計算科学をはじめよう！ ニュートンの冷却法則①
直腸温の連続測定による死後経過時間の推算（早大院先進理工）○（学）井上 幹康*（正）酒井 清孝（早大高等研）（正）山本 健一郎（防衛医大）（正）金武 潤

もし、あなたの目の前に亡骸があったら？

(引用;名探偵コナン「甘く冷たい宅配便」)
コナン君が遺体を見つけた時、かなり最初の段階で死後経過時間について考えますよね？死亡推定時刻が分かれば大きく容疑者を絞り込むことができますしね＼(^o^)／
そこで今回は死後経過時間について詳しい求め方を紹介します！

例題

室温16℃に設定された冷蔵車に男性の遺体があった。男性の体重は86kg、コナン君が遺体の直腸内温度を計ったところ、27℃であった。男性の死後経過時間を求めよ。

ニュートンの冷却の法則

まず、簡略化した事象を考えてみましょう。コップに入れられた熱いコーヒーが周囲の温度（空気）によって冷まされる問題を考えます。この時、次の「ニュートンの冷却の法則」が近似的に成り立つと言われています。
ニュートンの冷却の法則を表す微分方程式は

$\frac{dT}{dt} = -γ(T - T_s)$ (1)

で表されます。ここで、各文字は
$T$：湯の温度
$T_s$：周囲の温度
$\gamma$：冷却定数
を表しています。この微分法方程式は解析的に解くことができます。まず、(1)式の両辺を$T-T_s$で割ります。
$\frac{1}{T-T_s}dT = -γdt$

両辺を積分して積分定数をCとすると、
$log (T - T_s) = -γt + C$
となります。この式をTについて解くと下式のようになります。
$T = T_s + e^{-γt + C}$ (2)

初期条件$T(0)=T_0$を課すと、積分定数Cは
$C = log (T_0-T_s)$

となります。これを(2)式に代入すれば、解析解は以下で与えられます。
$T(t) = T_s + (T_0 - T_s) e^{-γt}$

（また、これを変形して）
$\frac{T(t)-T_s}{T_0-T_s} = e^{-γt}$

直腸温の連続測定による死後経過時間の推算

次に、ニュートンの冷却の法則を利用して死後経過時間を求めます。「直腸温の連続測定による死後経過時間の推算（早大院先進理工）」によると遺体の場合、微分方程式は次の様に変形されます。
$\frac{T_r-T_e}{T_o-T_e}=1.25e^{B*t}-0.25e^{5B*t}$ $(T_e\le23.2)$
$\frac{T_r-T_e}{T_o-T_e}=1.11e^{B*t}-0.11e^{10B*t}$ $(T_e\ge23.2)$
$B=-1.2815(bW)^{-0.625}+0.0284$
ただし、
$T_r$は直腸温℃
$T_e$は環境温℃
$T_o$は死亡時の直腸温であり37.2℃とする
$T_e$は環境温℃で測定時間の平均値で一定とする
$W$死体の体重kg
$b=1$は補正係数
$B$はNewtonの冷却定数
$t$は死後経過時間(h)
とします。
これを元に早速今回の例題を解きましょう！

import sympy

def death_time(Tr,Te,W):
    Tr=Tr
    Te=Te
    W=W
    B = sympy.Symbol('B')
    t = sympy.Symbol('t')
    To = 37.2
    b = 1
    E = sympy.S.Exp1
    equation1 = 1.25*(E**(B*t))-0.25*(E**(5*B*t))-(Tr-Te)/(To-Te)
    equation2 = 1.11*(E**(B*t))-0.11*(E**(10*B*t))-(Tr-Te)/(To-Te)
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    if Te <= 23.2:
        print(sympy.solve([equation1,equation3]))
    else:
        print(sympy.solve([equation2,equation3]))    
death_time(27,16,86) #今回の問題の設定       

$ $参考文献によると$T_e$は環境温の平均を使いますが、今回は発見時の室温の16℃を$T_e$としました。
結果は言わなくても分かると思いますが、scipyでこんな累乗変数を含む計算を直接できるわけもなく、朝起きても計算を続けていました(泣)。参考文献通りに試行錯誤法の考えの元で解くことにしました。

(引用Krylov Subspace Accelerated Newton Algorithm: Application to Dynamic Progressive Collapse Simulation of Frames)

関数の微分を求め方程式を反復の都度解く処理が大きい為、上図の様に微分を差分で近似させたクリロフ部分交換法を用いて計算することにしました。

import math
from scipy.optimize import newton_krylov
import sympy

#(3)式に相当する。体重(kg)からNewtonの冷却定数Bを算出。To(死亡時直腸温)は37.2とする
def predict_B(W):
    W=W
    B = sympy.Symbol('B')
    b = 1
    equation3 = -1.2815*((b*W)**(-0.625))+0.0284-B
    B = sympy.solve([equation3])[B]
    return B  

#(1)(2)式に相当。環境温Teで場合分け
def predict_deathtime(W,Tr,Te,ta):
    Bx = predict_B(W)
    Tr,Te,ta = Tr,Te,ta
    E = sympy.S.Exp1
    To = 37.2
    def F(x):
        if Te <= 23.2:
            return 1.25*(E**(Bx*x))-0.25*(E**(5*Bx*x))-(Tr-Te)/(To-Te)              
        else:
            return 1.11*(E**(Bx*x))-0.11*(E**(10*Bx*x))-(Tr-Te)/(To-Te)
    guess = ta
    sol = newton_krylov(F, guess, method='lgmres')
    print(sol)

#(体重,直腸温,環境温,死後経過時間の初期値)の順で入力
#死後経過時間の初期値を0によって微分値が求められないので注意する
predict_deathtime(86,27,16,2)

出力結果17.19191007875492

よって、今回の例題では17時間ほど前に死亡したことが分かりました！

参考文献

計算科学をはじめよう！ ニュートンの冷却法則①
直腸温の連続測定による死後経過時間の推算（早大院先進理工）○（学）井上 幹康*（正）酒井 清孝（早大高等研）（正）山本 健一郎（防衛医大）（正）金武 潤

概要

絶対パスから相対パスを取得するために以下のようなコードを実行した

import pathlib

p = pathlib.Path()
file_path = 'image-db'
file_path_rel = p.cwd().relative_to(file_path)

エラー

ValueError: '/directory/of/python' does not start with 'image-db'

relative_to()の中身がカレントディレクトリより外にあると動かないらしい．

対策

os.path.reipath('行先', '起点')を使えば複雑な相対パスも取得できる．

file_path_rel = os.path.relpath(file_path, os.getcwd())

参考サイト

https://hibiki-press.tech/python/os_path_abspath/1021

よくわかりませんが、

brew install libomp

pip uninstall lightgbm

pip install lightgbm

で解決しました。

PythonHackとは

本格的にやる前の準備段階

　動機づけ・・・Pythonこれからいいようだ。Web開発も本格的に使われていくと見込まれる。（と信じる）

コンテンツ

・Python　スタートブック　辻　真吾　著　技術評論社　を学ぶ
・PyQ　　　Python　はじめの一歩　オンラインレッスン　を学ぶ
・シラバス作成

シラバス

・変数
・計算・IF文
・list・for文
・辞書

初めに

• 本記事は、python・機械学習等々を独学している小生の備忘録的な記事になります。
• 「自身が気になったコードを写経しながら勉強していく」という、きわめてシンプルなものになります。
• 建設的なコメントを頂けますと幸いです（気に入ったらLGTM & ストックしてください）。

お題：Classify_images_Using_Python & Machine Learning

今日のお題は、Classify_images_Using_Python & Machine LearningというYoutube上の動画です。犬やら猫らの画像を学習させてそれを判定する、といったものです。

Classify Images Using Python & Machine Learning

分析はyoutubeの動画にある通り、Google Colaboratryを使用しました。

それではやっていきたいと思います。

　Step1: ライブラリのインポート~データの加工

import tensorflow as tf
from tensorflow import keras
from keras.models import Sequential
from keras.layers import Dense, Flatten, Conv2D, MaxPool2D, Dropout
from tensorflow.keras import layers
from keras.utils import to_categorical
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

次に

1. データの読み込み
2. データ型の確認
3. データシェイプ確認
4. いくつかのデータを確認
5. 画像を一つ出力

までをやっていきます。

#1
from keras.datasets import cifar10
(x_train, x_test), (y_train, y_test) = cifar10.load_data()

#2
print(type(x_train))
print(type(y_train))
print(type(x_test))
print(type(y_test))

#3
print('x_train shape:', x_train.shape)
print('y_train shape:', y_train.shape)
print('x_test shape:', x_test.shape)
print('y_test shape:', y_test.shape)

#4
index = 10
x_train[index]

#5
img = plt.imshow(x_train[index])

何やら画像が見えてきました。次にこの画像に紐づいているラベルを見ていきます。

print('The image label is:', x_train[index])

これを見ると4というラベルがついていることがわかります。次に開設されてますが、このデータセットは10種類の画像が含まれています。

# Get the image classification
classification = ['airplane', 'automobile', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck']
# Print the image class
print('The image class is:', classification[y_train[index][0]])

なので、これは「鹿」の画像だそうです（よく画面にめっちゃ寄っても厳しい...）。

次に被説明変数yをone_hot_encodingで、画像ラベルをNeural Networkに入れられる形に0 or 1(正しければ1、そうでなければ0)の数字を当てます。

y_train_one_hot = to_categorical(y_train)
y_test_one_hot = to_categorical(y_test)
print(y_train_one_hot)
print('The one hot label is:', y_train_one_hot[index])

次にデータセットを標準化(0 - 1)にします。

x_train = x_train / 255
x_test = x_test / 255
x_train[index]

これで凡そ加工はおしまいです。次にCNNモデルを組んでいきましょう！

Step2: モデルを構築

model = Sequential()

model.add( Conv2D(32, (5,5), activation='relu', input_shape=(32,32,3)))

model.add(MaxPool2D(pool_size=(2,2)))

model.add( Conv2D(32, (5,5), activation='relu'))

model.add(MaxPool2D(pool_size=(2,2)))

model.add(Flatten())

model.add(Dense(1000, activation='relu'))

model.add(Dropout(0.5))

model.add(Dense(500, activation='relu'))

model.add(Dropout(0.5))

model.add(Dense(250, activation='relu'))

model.add(Dense(10, activation='softmax'))

次にモデルをコンパイルし、学習させていきます。

model.compile(loss= 'categorical_crossentropy',
              optimizer= 'adam',
              metrics= ['accuracy'])


hist = model.fit(x_train, y_train_one_hot,
                 batch_size = 256,
                 epochs = 10,
                 validation_split = 0.2)

出来上がったmodelをtestデータでテストします。

model.evaluate(x_test, y_test_one_hot)[1]
>> 0.6811000108718872

うーん、動画でもそうでしたがあまり良い結果ではないですね。。。

とりあえずmatplotlibで精度と損失誤差を描画しておきます。

# Visualize the model accuracy
plt.plot(hist.history['accuracy'])
plt.plot(hist.history['val_accuracy'])
plt.title('Model accuracy')
plt.ylabel('Accuracy')
plt.xlabel('Epoch')
plt.legend(['Train', 'Val'], loc='upper left')
plt.show()

#Visualize the models loss
plt.plot(hist.history['loss'])
plt.plot(hist.history['val_loss'])
plt.title('Model loss')
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(['Train', 'Val'], loc='upper right')
plt.show()

最後に適当な画像（今回はネコの画像）を使ってモデルで予測してみます。

# Test the model with an example
from google.colab import files
uploaded = files.upload()

# show the image
new_image = plt.imread('cat-xxxxx.jpeg')
img = plt.imshow(new_image)

# Resize the image
from skimage.transform import resize
resized_image = resize(new_image, (32,32,3))
img = plt.imshow(resized_image)

# Get the models predictions
predictions = model.predict(np.array([resized_image]))
# Show the predictions
predictions

# Sort the predictions from least to greatest
list_index = [0,1,2,3,4,5,6,7,8,9]
x = predictions

for i in range(10):
  for j in range(10):
    if x[0][list_index[i]] > x[0][list_index[j]]:
      temp = list_index[i]
      list_index[i] = list_index[j]
      list_index[j] = temp

# Show the sorted labels in order
print(list_index)

# Print the first 5 predictions
for i in range(5):
  print(classification[list_index[i]], ':', round(predictions[0][list_index[i]] * 100, 2), '%')

結果は、

cat : 51.09 %
dog : 48.73 %
deer : 0.06 %
bird : 0.04 %
frog : 0.04 %

とネコかイヌかほぼ判定できていない結果に（笑）モデルの精度がイマイチ且つ使った画像が良くなかったです。

最後に

今回はtensorflowとkerasを使った画像判定の勉強をしました。もちろん使い物になる精度出ないことは重々承知ですが、足掛かりとしては良かったかなと思います。

引き続きよろしくです。

(これまでの学習)
1. 【1日1写経】Predict employee attrition【Daily_Coding_001】
2. 【1日1写経】Build a Stock Prediction Program【Daily_Coding_002】

はじめに

一度は勉強したけど、「あれ？これってどうだったっけ...」という方向けのまとめです。

FE（基本情報技術者試験）でもPythonが選択できるようになりました。
本来なら今年の春試験が初の出題となる予定でしたがコロナの影響で中止となりましたね。

そんなPythonのlistもFEの例題でバッチリ出題されているので、使い方をまとめてみました。

試験の直前に見返すも良し、普段のプログラミングで忘れがちなメソッドを思い出すも良し。
ぜひご活用ください。

listの初期化

まずは基本のlistの宣言から。

alist = []
blist = [1, 1, 2, 3, 5, 8, 13]
clist = list(range(5))
dlist = [i for i in range(10) if (i % 2 == 0)]  # リスト内包表記

print(f"alist = {alist}\nblist = {blist}\nclist = {clist}\ndlist = {dlist}")

output

alist = []
blist = [1, 1, 2, 3, 5, 8, 13]
clist = [0, 1, 2, 3, 4]
dlist = [0, 2, 4, 6, 8]

代入とコピー

意外と盲点です。なので、以降のlistの扱いでは参照渡しと値渡しの両方を交えて説明します。

• = でlistを代入
  • 参照渡し
  • 代入元(alist)が変更 → blist も変わる
• .copy() メソッドの返り値を代入
  • 値渡し
  • 代入元(alist)が変更 → clist は変化無し

alist = [1, 2, 3]

# 参照渡し
blist = alist

# 値渡し
clist = alist.copy()

def delete_middle(list_obj, msg):
    list_obj.pop(int(len(list_obj)/2))
    print(msg, list_obj)

delete_middle(alist, "alist =")
print(f"blist = {blist}")
print(f"clist = {clist}")

output

alist = [1, 3]
blist = [1, 3]
clist = [1, 2, 3]

listに要素を追加

• 要素を後ろに1つだけ追加
  • .append(x) メソッド
• 指定したindexの場所に1つ追加
  • .insert(index, obj) メソッド
• listにlistを追加
  • + operator(演算子)を使う
  • 演算後に blist や clist の要素が変化しても dlist には影響なし

alist = [1, 2, 3]
blist = alist
clist = alist.copy()

# 4を後ろに追加
alist.append(4)

# 1のindexに100を追加
alist.insert(1, 100)

# listにlistを追加
dlist = blist + clist

print(f"alist = {alist}\nblist = {blist}\nclist = {clist}\ndlist = {dlist}")

alist = [1, 100, 2, 3, 4]
blist = [1, 100, 2, 3, 4]
clist = [1, 2, 3]
dlist = [1, 100, 2, 3, 4, 1, 2, 3]

listの中身を削除

indexを指定して削除

alist = [1, 3, 5, 7, 9]
blist = alist
clist = alist.copy()

del alist[2]
print(f"del alist[2] =>\n\talist = {alist}\n\tblist = {blist}\n\tclist = {clist}\n")

alist.pop(2)
print(f"alist.pop(2) =>\n\talist = {alist}\n\tblist = {blist}\n\tclist = {clist}")

output

del alist[2] =>
    alist = [1, 3, 7, 9]
    blist = [1, 3, 7, 9]
    clist = [1, 3, 5, 7, 9]

alist.pop(2) =>
    alist = [1, 3, 9]
    blist = [1, 3, 9]
    clist = [1, 3, 5, 7, 9]

値を指定して削除

先頭側にある値だけが削除されます。その値全てではありません。

alist = [1, 2, 2, 3, 3, 3]
blist = alist
clist = alist.copy()

alist.remove(2)

print(f"alist = {alist}\nblist = {blist}\nclist = {clist}")

output

alist = [1, 2, 3, 3, 3]
blist = [1, 2, 3, 3, 3]
clist = [1, 2, 2, 3, 3, 3]

一番後ろの項目を削除

.pop() に引数を入れなければ .pop(-1) と同じで、リストの最後の要素が消されます。

alist = [1, 3, 5, 7, 9]
blist = alist
clist = alist.copy()

alist.pop()
print(f"alist.pop(2) =>\n\talist = {alist}\n\tblist = {blist}\n\tclist = {clist}")

output

alist.pop() =>
    alist = [1, 3, 5, 7]
    blist = [1, 3, 5, 7]
    clist = [1, 3, 5, 7, 9]

全て削除

alist = [1, 2, 3]
blist = alist
clist = alist.copy()

# alist = [] もやってることは同じ。
# alist = [] だと宣言と紛らわしいけど、開発チームのお好みで。
alist.clear()

print(f"alist = {alist}\nblist = {blist}\nclist = {clist}")

output

alist = []
blist = []
clist = [1, 2, 3]

値の順番で並び替える

alist = [20, 1, 5, 13, 8]
blist = alist
clist = alist.copy()

alist.sort()

print(f"alist = {alist}\nblist = {blist}\nclist = {clist}")

output

alist = [1, 5, 8, 13, 20]
blist = [1, 5, 8, 13, 20]
clist = [20, 1, 5, 13, 8]

最初に現れる値のindexを調べる

alist = [333, 1, 333, 22, 333, 22]
# 22が何番目にあるか
index = alist.index(22)
print(f"index = {index}")

output

index = 3

値が何個含まれているか

alist = [0, 1, 1, 0, 1, 0 ,1, 0, 0]
# 1が何回出てきたか
count = alist.count(1)
print(f"count = {count}")

output

count = 4

さいごに

メジャーどころのlistのメソッドは抑えたかなと思います。
参照渡しと値渡しの挙動を説明したくて返って読みにくくなってしまいました。すみませんmm

忘備録としてお使いいただければと幸いです。それでは。

結論

MacにPythonをインストールして、Pythonを使えるようにする。
MacはPython２系が元々インストールされているが今回はPython３系をインストールする。

背景

Pythonを始めるにはまずはPythonを使えるようにしなければならない。プログラミングを始めるにあたって一番はじめにつまずくのが環境構築。これからPythonを始めたいと思っている人の参考になるように手順を書きたいと思ったため。初心者向けのため仮想環境にインストールするなど複雑なことは考えない。

全体の流れ

① コマンド操作を覚える
② Pythonをインストールする

コマンド操作を覚える

Pythonをインストールする前にMacのコマンド操作を覚えること。
これが少しできるだけでプログラミングも楽しくなる。

ls　　　　# フォルダ内を表示する
ls -a　　# 隠しファイルも含めファイルやフォルダを全て表示する
cd　　　　# フォルダを移動する
pwd　　　# フォルダのパスを表示する
touch　　# 空のファイルを作る
mkdir　　# フォルダを作る
mv　　　　# フォルダ・ファイルを移動する
cp　　　　# コピーする
rm　　　　# ファイルを削除する
rm -r　　# フォルダを削除する
open　　　# ファイルを開く

とりあえずこれだけ覚えておけばなんとなく操作できると思う。

Pythonをインストールする

1. 公式サイトからPythonをダウンロードする。

https://www.python.org/downloads/

無事にインストールできたら下記のコマンドをターミナルに入力する。
インストールしたPythonのバージョンが表示されればOK。

python3 --version

1. PATHを設定する（PATHを通す）

.bash_profile（ホームディレクトリの直下に隠しファイルになっている）を確認する。
下記のような内容が記載されているはず。なければファイルを作成し、同じように書く。（PATHを通す）

# Setting PATH for Python 3.8
# The original version is saved in .bash_profile.pysave
PATH="/Library/Frameworks/Python.framework/Versions/3.8/bin:${PATH}"
export PATH

.bashrcに下記を記載する。なければファイルを作成し、同じように書く。

alias python=python3
alias pip=pip3

さらに、.bash_profileに下記を追記する。

#.bashrc読み込み用
source ~/.bashrc

これでPATHを通すことができた。
MacではPython２系がインストールされているので、PATHを通すと今後の開発がスムーズになる。
いちいちPython３を入力しなくてもPythonで実行できるようになる。

まとめ

プログラミングではじめにつまずくのが環境構築だ。PATHを通すっていうところでつまずくと思う。今回は細かい内容は書いていないが、とにかく手順どおりにやってみてほしい。これができれば、Pythonを始められる。

参考にした情報

• Macのターミナルコマンド一覧（基本編）
• 【Python3】まっさらなMacにPythonの最新バージョンをインストールする手順まとめ

はじめに

　私は、時々Flaskを使うのですが、使い方をよく忘れてしまうので、ここにメモをしておきます。

レスポンスを受け取るまで

インストール

　まず、Flaskをインストールします。以下のコマンドを打ち込むだけで良いです。

$ pip install Flask

ファイル

　次に適当にサーバの動作を行うファイルを作ります。以下に示すのは、必要最低限の動作を行うものです。

app.py

from flask import Flask

app = Flask(__name__)

@app.route("/")
def hello():
    return "hello"

if __name__ == "__main__":
    app.run(debug=True)

サーバの稼働

　そして次のコマンドをこのファイルが存在するディレクトリ上で打ち込みます。

$ python app.py

これでhttp://localhost:5000にアクセスするとhelloと帰ってくるようになりました。

基本の関数の使い方

@app.route(URL, methods)

これはその後に来る関数をリクエストが来た際に走らせるようにするラッパーです。リクエストを受け取る場合は必ず使います。引数はそれぞれ、
URL: リクエストを受け取るURLのroute
methods: リクエストのメソッド
使い方としては以下のようなものがあります。また、URLから直接変数を得る方法などもあります。

基本的な使い方

# http://localhost:5000/hogeにGETリクエストを送ることで動く
@app.route("/hoge", methods=["GET"])
def hoge():
    # ここで処理が行われる

URLから変数を取得する

FlaskではURLの文字列をデータとして取り込むことができます。例えば以下の例では、http://localhost:5000/hoge/fugaに対してGETメソッドを送ると、「hello! fuga」と返ってきます。

@app.route("/hoge/<text>", methods=["GET"])
def hello(text):
    return "hello! " + text

app.run(host, path, debug)

　ここではアプリの稼働を行います。それぞれの引数は、
host: アプリを稼働するipアドレス、デフォルトはlocalhost
port: アプリにアクセスできるポート、デフォルトは5000
debug: Trueの時はアプリ稼働時にデバッグ出力を行う

これらを使って、変更を加えるとしたら、以下のような感じになります。

app.run(host="0.0.0.0", port=8080, debug=True)

hostを上のように設定しておくことで、外部のデバイスからでもこのコードを走らせているPCのアドレスが分かれば、アクセスできるようになります。

jsonの取り扱い

リクエストで受け取る方法

　リクエストでjsonを受け取る時には、jsonモジュールを使用します。

from flask import Flask
import json

@app.route("/", methods=["POST"])
def hoge():
    req_data = json.loads(request.data)
    return req_data["hoge"]

json.loads関数を用いることで、jsonのデータを配列や辞書などのpythonで扱える形に変換しています。

レスポンスで出力する方法

　レスポンスにjsonを含める場合には、Flaskのjsonify関数を使用します。

from flask import Flask, jsonify

@app.route("/", methods=["GET"])
def hoge():
    res_data = {"json_data": "hogehoge"}
    return jsonify(res_data)

　基本的にはデータをjsonifyで包んで送るだけです。

最後に

　何か間違いや、足りない部分などがあれば教えていただけるとありがたいです。また、質問なども受け付けます。　

はじめに

これは自分用の備忘録なのでご了承を。

Nanについて

Nanは不動点少数における考え方であり、異常な数値であることを表している。
欠損値という名前がついておるようだ。
イメージは、エクセルで式を計算したときにでるエラーでいいかな。
Noneではないことに注意する。

Nanの性質

Nanには重要な性質がある。
性質には、以下のものがある。

・Nanを含む四則演算の結果はNanが返ってくる
・Nanの比較演算子は必ずFalseとなる

通常の計算とは結果が異なるため、注意が必要。

numpyのisnan()関数

numpyに入っているisnan()はNanがあるかどうかを判定する関数。
そのため、if文などでこれを使うと、Trueが返ってくる。

nansum()関数

欠損値を含む計算はNanが返ってきてしまうのだ。
けど、Nanを含む数値群を計算したい！という場合はnansum()を使う。
これを使うことで、Nanを除いた値の合計が計算される。
便利ねーーー。

nan_to_num()関数

Nanを別のものに変えたいーーって時のための関数。
デフォルトでは、0に変更されて元のndarrayは変更されないみたい。
第一引数で変更する配列を指定して、第二引数で元の配列を保持するか指定、第三引数で置換する文字、数字を指定。

はじめに

AtCoder Problems の Recommendation を利用して、過去の問題を解いています。
AtCoder さん、AtCoder Problems さん、ありがとうございます。

今回のお題

AtCoder Regular Contest D - Grid Coloring
Difficulty: 855

今回のテーマ、シミュレーション

内容は難しいところはないと思いますが、実装に手間がかかりそうな問題です。

1 2 2 3 3
4 4 4 4 3
5 5 5 5 5

与えられた数値で上から順に左右に塗りつぶせばOKです。

Ruby Array

ruby.rb

h, w = gets.split.map(&:to_i)
_ = gets.to_i
a = gets.split.map(&:to_i)
m = Array.new(h){Array.new(w, 0)}
y, x, lr = 0, 0, 1
a.each_with_index do |n, v|
  n.times do
    m[y][x] = v + 1
    x += lr
    if x == w
      y += 1
      x = w - 1
      lr = -1
    elsif x < 0
      y += 1
      x = 0
      lr = 1
    end
  end
end
h.times do |i|
  puts m[i].join(' ')
end

join.rb

  puts m[i].join(' ')

  puts m[i] * ' '

最近知ったのですが、joinは上記の様に書くこともできます。

Ruby Class

ruby.rb

class MASS
  def initialize(h, w)
    @h = h
    @w = w
    @m = Array.new(@h){Array.new(@w, 0)}
    @x = 0
    @y = 0
    @LR = 1
  end
  def draw(v, n)
    n.times do
      @m[@y][@x] = v + 1
      @x += @LR
      if @x == @w
        @y += 1
        @x = @w - 1
        @LR = -1
      elsif @x < 0
        @y += 1
        @x = 0
        @LR = 1
      end
    end
  end
  def out
    @h.times do |i|
      puts @m[i].join(' ')
    end
  end
end

h, w = gets.split.map(&:to_i)
_ = gets.to_i
a = gets.split.map(&:to_i)
m = MASS.new(h, w)
a.each_with_index do |v, i|
  m.draw(i, v)
end
m.out

今回の問題では、あまり有効ではないのですが、学習を兼ねてclassを使用してみました。
コード長が倍近くになりましたが、実行時間は変わらないようです。

Python

python.py

from sys import stdin

def main():
    input = stdin.readline
    h, w = map(int, input().split())
    _ = int(input())
    a = list(map(int, input().split()))
    m = [[0] * w for _ in range(h)]
    x, y, LR, v = 0, 0, 1, 0
    for n in a:
        v += 1
        for _ in range(n):
            m[y][x] = str(v)
            x += LR
            if x == w:
                y += 1
                x = w - 1
                LR = -1
            elif x < 0:
                y += 1
                x = 0
                LR = 1
    for i in range(h):
        print(" ".join(m[i]))
main()

Python のコード長が長いのは、半角スペースもカウントしているからでしょうか。

	Ruby Array	Ruby Class	Python
コード長 (Byte)	392	631	603
実行時間 (ms)	17	18	25
メモリ (KB)	2428	3836	3828

まとめ

• ARC 080 D を解いた
• Ruby に詳しくなった
• Python に詳しくなった

パラメータの最適化

実測データとモデルとの比較を行うとき、モデルの関数に含まれるパラメータの最適化は必須です。

大学に通っていた頃、自分の手でパラメータをちょっとずつ変えて最適化をしている努力家がいました。ソフトウェアが専門でない学科においては多少仕方のない光景ですが、やはりこういうのはプログラミングで解決しましょう。

今回紹介する2手法

• random search method (ランダムサーチ法)
• gradient descent method (最急降下法)

そもそもどのように最適化するのか

最適化プログラムの構成

これを何万回と繰り返し、最適なパラメータへと更新し続けます。

評価の仕方

最も一般的なのは、最小自乗法です。
今回の私のプログラムにも採用しています。

誤差をそのままではなく2乗してから足すことで、実測とモデルの誤差を積算することができます。

random search method (ランダムサーチ法)

パラメータにランダムな値を入れて，評価関数が一番小さくなる値を見つける手法です。

gradient descent method (最急降下法)

評価関数の偏微分値を求め、評価関数が小さくなる向きに少しずつパラメータをずらしていく手法です。

　2手法の比較

解の精度は、ランダムサーチ法　<<　最急降下法

かなり精度に差があります。
ランダムサーチでは最適解に近づけば近づくほど、より最適なパラメータに近づくことが難しくなります。それに対して、最急降下法では堅実に最適解に近づくことができます。

ランダムサーチ法の試行回数をかなり増やしたとしても、再急降下法より良い精度を出すことはほぼ不可能でしょう。

参考程度に、4パラメータの最適化で10万回の試行回数でプログラムを回した結果を見てみましょう。

悲しいことに、10万回ランダムサーチしても、毎回答えが違うしグラフは全然沿いません。
精度面は、最急降下法の圧勝ですね。

最急降下法の欠点

そんな最急降下法にも致命的な欠点があります。

局所的な最適解にはまってしまうと、さらなる最適解があっても見つけられなくなるのです。

コンピュータに全てを委ねる最適化プログラムにおいて、「最適じゃない解」を出力する可能性があるのは見過ごせません。

2手法のコンビネーション

そこで、最も良い解決法が、「ランダムサーチ法と最急降下法のコンビネーション」です。

1. 大きな範囲を指定してランダムサーチをする。
2. 繰り返し行うことで最適な値があると推測される範囲を削っていく。
3. 最終的に、右に記した範囲になるまでランダムサーチを行う。
4. 最急降下法により最適解を求める。

これで問題解決です。

最後に、2手法のコード

実測値データは私が回路の周波数応答の計測をしたときのものです。

ランダムサーチ法

random_search.py

# random search method

import random
import numpy as np
import math
import matplotlib.pyplot as plt


#----- data -----

freq = np.array([0.05, 0.1, 0.2, 0.4, 0.6061, 0.8, 1, 10,100])

ReR = np.array([0.005890792, 0.005824175, 0.005726238, 0.005501362, 0.005399829, 0.005244864, 0.0051687, 0.004736195, 0.004552081])

ImR = np.array([0.00476029, 0.002511934, 0.001470715, 0.000984505, 0.000808862, 0.000718529, 0.000655644, 0.00021882, -0.0000045569]) 

num=100000 # number of trials

#----------------




lsv = 100000
Flist=[0 for i in range(num)]
R1list=[0 for i in range(num)]
R2list=[0 for i in range(num)]
C1list=[0 for i in range(num)]
C2list=[0 for i in range(num)]
N = len(freq)
R_opt = np.array([0j for i in range(N)])
theo = np.array([0j for i in range(N)])
R1_opt,R2_opt,C1_opt,C2_opt = 0,0,0,0



def init():
    sv = 0 # squares value
    R1 = random.uniform(0.001, 0.01)
    R2 = random.uniform(0.001, 0.01) 
    C1 = random.uniform(0, 500)
    C2 = random.uniform(0, 500)
    return sv,R1,R2,C1,C2


def func_theoretical(R1,R2,C1,C2,freq):
    c1 = 1/(2j * math.pi * freq * C1)
    c2 = 1/(2j * math.pi * freq * C2)
    c3 = R2 + c2
    theo = R1 + c1*c3/(c1+c3)
    return theo


def func_distance(ReR, ImR, theo):
    return (ReR - theo.real) **2 + (ImR + theo.imag) **2


#--- traials ---
for i in range(num):
    sv,R1,R2,C1,C2 = init()

    for i in range(N):
        theo[i] = func_theoretical(R1,R2,C1,C2,freq[i])
        sv += func_distance(ReR[i], ImR[i], theo[i])

    if sv < lsv:
        lsv = sv
        R_opt = theo
        R1_opt,R2_opt,C1_opt,C2_opt = R1,R2,C1,C2
#----------------


#--- result ---
print('number of trials :',num)
print('Least Squares Value :',lsv)
print('R1,R2,C1,C2 :' ,R1_opt,R2_opt,C1_opt,C2_opt)


#--- graph ---
plt.title('random search', loc='center')
plt.plot(ReR,ImR, label="Measured value")
plt.plot(R_opt.real, - R_opt.imag ,label="Theoretical value")
plt.grid(color='gray')
plt.legend(loc="upper left")
plt.show()

最急降下法

gradient_descent.py

import random
import numpy as np
import math
import matplotlib.pyplot as plt


#--- data ---

freq = np.array([0.05, 0.1, 0.2, 0.4, 0.6061, 0.8, 1, 10,100])

ReR = np.array([0.005890792, 0.005824175, 0.005726238, 0.005501362, 0.005399829, 0.005244864, 0.0051687, 0.004736195, 0.004552081])

ImR = np.array([0.00476029, 0.002511934, 0.001470715, 0.000984505, 0.000808862, 0.000718529, 0.000655644, 0.00021882, -0.0000045569]) 

#------------


#--- setting ---

ParaList = [0.005, 0.0014, 138 , 546] #R1,R2,C1,C2

d = [0.00001, 0.000001, 0.1, 0.1] #d(R1,R2,C1,C2)

num=1000 # number of trials
pf=11 # print freqency

#----------------

def func_lsv():
    global freq, ReR, ImR, ParaList

    R1,R2,C1,C2 = ParaList
    for Freq, rer, imr in zip(freq, ReR, ImR):
        a = 1/(2j * math.pi * freq * C1)
        b = 1/(2j * math.pi * freq * C2)
        c = R2 + b
        R = R1 + a*c/(a+c) 
        distance = (ReR - R.real) **2  + (ImR + R.imag) **2
    lsv = sum(distance)
    return lsv


def gradient_descent(N,Rs):
    global ParaList

    ParaList[N] = ParaList[N] - Rs
    lsv_under = func_lsv()
    ParaList[N] = ParaList[N] + 2 * Rs
    lsv_over = func_lsv()

    if lsv_under < lsv_over :
        ParaList[N] = ParaList[N] - 2 * Rs
        lsv = lsv_under
    else:
        lsv = lsv_over

    return ParaList, lsv


#--- main ---
cnt_pf=0
print("  R1  /  R2  /  C1  /  C2  /  Least Squares Value")

for i in range(num):
    cnt_pf +=1
    for N in range(4):
        ParaList, lsv = gradient_descent(N, d[N])
    if cnt_pf == pf:
        print(ParaList,lsv)
        cnt_pf=0


Rbest=np.array([0j for i in range(len(freq))])

R1,R2,C1,C2 = ParaList
for i in range(len(freq)):
    a = 1/(2j * math.pi * freq[i] * C1)
    b = 1/(2j * math.pi * freq[i] * C2)
    c = R2 + b
    Rbest[i] = R1 + a*c/(a+c) 


#--- result ---
print('trials :',num)
print('Least Squares Value: :', lsv)
print('R1 ,R2, C1, C2 :',ParaList)


#--- graph ---
plt.title('saikyuukouka', loc='center')
plt.plot(ReR,ImR, label="Measured value")
plt.plot(Rbest.real, -Rbest.imag ,label="Theoretical value")
plt.grid(color='gray')
plt.legend(loc="upper left")
plt.show()