物理の駅 by 現役研究者

Physics Station → PhSt 質問・疑問・間違いの指摘は、コメントに書くか、直接伝えるときっと良いことがあります。主にWindows or Ubuntu用の記事です

Python で数値を科学的表記にする方法

科学的な表記のため、あまりにも長い桁数を丸める必要がある。ざっくりそういう需要のための関数を作った。

やり方は、指数表記にして指数部の大きさで判定させているだけである。自分で作ったほうが早いと思うので参考にしつつ作ってみてはどうだろうか。

def align_digit(f, digit=None,shift=3):
    if digit==None:
        digit = int("{:e}".format(f).split("e")[1])-shift
    if digit<=0: s = str("{:."+str(-digit)+"f}").format(f)
    else: s = str("{:.3e}").format(f)
    return s, digit

print(align_digit(0.00001))
print(align_digit(0.0000100123))

print(align_digit(0.0001))
print(align_digit(0.000100123))

print(align_digit(0.001))
print(align_digit(0.00100123))

print(align_digit(0.01))
print(align_digit(0.0100123))

print(align_digit(0.1))
print(align_digit(0.100123))

print(align_digit(1))
print(align_digit(1.00123))

print(align_digit(10))
print(align_digit(10.0123))

print(align_digit(100))
print(align_digit(100.123))

print(align_digit(1000))
print(align_digit(1001.23))

print(align_digit(10000))
print(align_digit(10012.3))

print(align_digit(100000))
print(align_digit(100123))

shift=3 桁分余分に表示するという感じ。結果は以下の通り

('0.00001000', -8)
('0.00001001', -8)
('0.0001000', -7)
('0.0001001', -7)
('0.001000', -6)
('0.001001', -6)
('0.01000', -5)
('0.01001', -5)
('0.1000', -4)
('0.1001', -4)
('1.000', -3)
('1.001', -3)
('10.00', -2)
('10.01', -2)
('100.0', -1)
('100.1', -1)
('1000', 0)
('1001', 0)
('1.000e+04', 1)
('1.001e+04', 1)
('1.000e+05', 2)
('1.001e+05', 2)

戻り値のtupleの2つめを再利用してdigit=に入れると、同じ桁数の表記が得られる。

Python scipyで任意の方程式の解を求める

from scipy.optimize import fsolve
import numpy as np

# 変数
beta = 0.6955076793404303
ionpair = 4886

# 方程式
func = lambda delE : np.sqrt(delE/(np.log(ionpair*beta**2)-np.log(1-beta**2)-beta**2))*beta-92

# 初期値
initial_guess = 1000

# 解く
solution = fsolve(func, initial_guess)

# 解
print(solution[0])

139021.52385380902 と出てくるはず。確認のためにグラフを作っておく

x = np.linspace(0,200000)
y = func(x)
import matplotlib.pyplot as plt
plt.plot(x, y)
plt.grid(True)
plt.xlabel("x")
plt.ylabel("y")
plt.show()

f:id:onsanai:20210709091435p:plain

Python matplotlibのhist2dでカラーバーをfigsなしで表示する方法

import matplotlib.pyplot as plt
import numpy as np

x = np.random.normal(0,1,10000) #平均0、標準偏差1、N=10000
y = np.random.normal(3,2,10000) #平均3、標準偏差2、N=10000

counts, xedges, yedges, im = plt.hist2d(x, y, range=[[-10,10],[-10,10]], bins=[100,100], cmap="binary")

plt.xlabel("X")
plt.ylabel("Y")
im.set_clim(0,50) #カラーバーの表示範囲
plt.colorbar(im,label="Entries") #カラーバーを表示しラベルをセット
plt.savefig("hist2d.pdf")
plt.show()

f:id:onsanai:20210707092815p:plain

カラーバーのスケールを対数(ログ)にする場合は、normを以下のように指定する。表示する範囲の下限値に0が使えないことに留意

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

x = np.random.normal(0,1,10000)
y = np.random.normal(3,2,10000)
counts, xedges, yedges, im = plt.hist2d(x, y, range=[[-10,10],[-10,10]], bins=[100,100], cmap="binary", norm=matplotlib.colors.LogNorm())

plt.xlabel("X")
plt.ylabel("Y")
im.set_clim(0.1,100)
plt.colorbar(im,label="Entries")
plt.savefig("hist2d_log.pdf")
plt.show()

f:id:onsanai:20210707093135p:plain

Pythonで、正規分布でフィッティングする汎用的な関数を作ってみる 2

過去のコードをより汎用的にするために修正した。

phst.hateblo.jp

小数の桁を取得して、小数点以下の有効桁数が4桁になるようにした。

gaus_sample.txtガウス分布のフィッティング用のサンプルデータ Sample data からダウンロードできます。

import math

def gaussian_func(x, constant, mean, sigma):
    import numpy as np
    return constant * np.exp(- (x - mean) ** 2 / (2 * sigma ** 2))

def align_digit(f, digit=None,shift=3):
    if digit==None:
        digit = int("{:e}".format(f).split("e")[1])-shift
    if digit<=0: s = str("{:."+str(-digit)+"f}").format(f)
    else: s = str(f)
    return s, digit

def fit_gaussian_impl(arr_x, arr_y, mean_value, stdev_value, show_plot, show_result):

    '''
    return chi2, constant, error, mean, error, sigma, error, degree of freedom, p-value
    '''

    arr_yerror = [math.sqrt(y) for y in arr_y]

    import numpy as np
    parameter_initial = np.array([max(arr_y), mean_value, stdev_value])
    
    arr_x2 = []
    arr_y2 = []
    arr_yerror2 = []
    for x, y, yerror in zip(arr_x, arr_y, arr_yerror): #entry=0のデータを除外
        if y == 0:continue
        arr_x2.append(x)
        arr_y2.append(y)
        arr_yerror2.append(yerror)

    if len(arr_x2) < 4:return {}

    from scipy.optimize import curve_fit
    popt, pcov = curve_fit(gaussian_func, arr_x2, arr_y2, sigma=arr_yerror2, absolute_sigma =True, p0=parameter_initial)
    stderr = np.sqrt(np.diag(pcov)) #対角行列を取って平方根

    arr_fitted_y = gaussian_func(np.array(arr_x2), popt[0], popt[1], popt[2]) 

    from scipy.stats import chisquare
    chisq, p = chisquare(f_exp=arr_y2, f_obs=arr_fitted_y, ddof = 2)
    ndf = len(arr_x2) - 3

    mat = np.vstack((popt,stderr)).T
    if show_result:
        import pandas as pd
        df = pd.DataFrame(mat,index=("Constant", "Mean", "Sigma"), columns=("Estimate", "Std. error"))
        print(df)

    if show_plot:
        import matplotlib.pyplot as plt
        plt.errorbar(arr_x,arr_y,linestyle="None",marker="+",yerr=arr_yerror,label="Data",color="tab:blue")

        arr_curve = gaussian_func(np.array(arr_x),popt[0],popt[1],popt[2])
        str_popt0, digit_popt0 = align_digit(popt[0])
        str_popt1, digit_popt1 = align_digit(popt[1])
        str_popt2, digit_popt2 = align_digit(popt[2])
        str_serr0, _ = align_digit(stderr[0], digit_popt0)
        str_serr1, _ = align_digit(stderr[1], digit_popt1)
        str_serr2, _ = align_digit(stderr[2], digit_popt2)
        plt.plot(arr_x,arr_curve,color="tab:orange",
                label=str("Fitted\n$\\chi^2$/ndf: {0:.2f}/{1}\n"+
                          "Constant: {2}$\\pm${3}\n"+
                          "Mean: {4}$\\pm${5}\n"+
                          "Sigma: {6}$\\pm${7}").format(chisq,ndf,
                                                    str_popt0,str_serr0,
                                                    str_popt1,str_serr1,
                                                    str_popt2,str_serr2))
        plt.legend(bbox_to_anchor=(1.12, 1.15), loc='upper right', borderaxespad=0)
        plt.savefig("fit.pdf")
        plt.show()
    
    obj = {}
    obj["chi2"]    =chisq
    obj["constant"]=[popt[0],stderr[0]]
    obj["mean"]    =[popt[1],stderr[1]]
    obj["sigma"]   =[popt[2],stderr[2]]
    obj["ndf"]     =ndf
    obj["pvalue"]  =p
    return obj

def fit_gaussian(vx, nbin, min_x, max_x,*,mean_value=None,stdev_value=None,show_plot=False,show_result=True):
    '''
    return chi2, constant, error, mean, error, sigma, error, degree of freedom, p-value
    '''
    if len(vx) == 0:return {}

    arr_x = [0 for i in range(nbin)]
    arr_y = [0 for i in range(nbin)]

    wbin = (max_x - min_x) / nbin
    if wbin <= 0:return {}

    for i in range(nbin): arr_x[i] = min_x + wbin * (i + 0.5)

    for x in vx:
        bin = math.floor((x - min_x) / wbin)
        if bin < 0:continue
        if bin >= nbin:continue
        arr_y[bin]+=1

    if mean_value == None:
        from statistics import mean
        mean_value = mean(vx)
    if stdev_value == None:
        from statistics import stdev
        stdev_value = stdev(vx)

    if show_plot:
        import matplotlib.pyplot as plt
        arr_yerror = [math.sqrt(y) for y in arr_y]
        str_mean_value, _ = align_digit(mean_value)
        str_stdev_value, _ = align_digit(stdev_value)
        plt.errorbar(arr_x,arr_y,linestyle="None",marker="+",yerr=arr_yerror,color="tab:blue",
                    label="Data\nEntries: {0}\nMean: {1}\nStdev: {2}".format(len(vx),str_mean_value,str_stdev_value))
        print(mean_value,stdev_value)
        plt.legend(bbox_to_anchor=(1.12, 1.15), loc='upper right', borderaxespad=0)
        plt.show()

        
    return fit_gaussian_impl(arr_x, arr_y, mean_value, stdev_value, show_plot, show_result)

import numpy.random
fit_gaussian(numpy.random.randn(100000),50,-5,5,show_plot=True)

arr_x=[]
arr_y=[]
with open("gaus_sample.txt") as f:
    lines = f.readlines()
    for line in lines:
        arr_x.append(float(line.split()[0]))
        arr_y.append(float(line.split()[1]))

fit_gaussian_impl(arr_x, arr_y, mean_value=0, stdev_value=1, show_plot=True, show_result=True)

最後に表示されるグラフ (詳細は Python3でROOTと同様にフィッティングとパラメータの標準誤差を算出する - 物理の駅 by 現役研究者 参照)

f:id:onsanai:20210707154716p:plain

C#でGoogleのワンタイムパスワード計算アルゴリズム TOTP (RFC6238)を実装する

Stackoverflowのこの投稿で十分

stackoverflow.com

だが、クリップボードにコピーして即プログラムを終了するように少し修正した。

using System;
using System.Windows.Forms;
using System.Security.Cryptography;

namespace totp
{
    class Program
    {
        public class Totp
        {
            const long unixEpochTicks = 621355968000000000L;
            const long ticksToSeconds = 10000000L;
            private const int step = 30;
            private const int totpSize = 6;
            private byte[] key;

            public Totp(byte[] secretKey)
            {
                key = secretKey;
            }

            public string ComputeTotp()
            {
                var window = CalculateTimeStepFromTimestamp(DateTime.UtcNow);

                var data = GetBigEndianBytes(window);

                var hmac = new HMACSHA1();
                hmac.Key = key;
                var hmacComputedHash = hmac.ComputeHash(data);

                int offset = hmacComputedHash[hmacComputedHash.Length - 1] & 0x0F;
                var otp = (hmacComputedHash[offset] & 0x7f) << 24
                       | (hmacComputedHash[offset + 1] & 0xff) << 16
                       | (hmacComputedHash[offset + 2] & 0xff) << 8
                       | (hmacComputedHash[offset + 3] & 0xff) % 1000000;

                var result = Digits(otp, totpSize);

                return result;
            }

            public int RemainingSeconds()
            {
                return step - (int)(((DateTime.UtcNow.Ticks - unixEpochTicks) / ticksToSeconds) % step);
            }

            private byte[] GetBigEndianBytes(long input)
            {
                // Since .net uses little endian numbers, we need to reverse the byte order to get big endian.
                var data = BitConverter.GetBytes(input);
                Array.Reverse(data);
                return data;
            }

            private long CalculateTimeStepFromTimestamp(DateTime timestamp)
            {
                var unixTimestamp = (timestamp.Ticks - unixEpochTicks) / ticksToSeconds;
                var window = unixTimestamp / (long)step;
                return window;
            }

            private string Digits(long input, int digitCount)
            {
                var truncatedValue = ((int)input % (int)Math.Pow(10, digitCount));
                return truncatedValue.ToString().PadLeft(digitCount, '0');
            }

        }
        public static class Base32Encoding
        {
            public static byte[] ToBytes(string input)
            {
                if (string.IsNullOrEmpty(input))
                {
                    throw new ArgumentNullException("input");
                }

                input = input.TrimEnd('='); //remove padding characters
                int byteCount = input.Length * 5 / 8; //this must be TRUNCATED
                byte[] returnArray = new byte[byteCount];

                byte curByte = 0, bitsRemaining = 8;
                int mask = 0, arrayIndex = 0;

                foreach (char c in input)
                {
                    int cValue = CharToValue(c);

                    if (bitsRemaining > 5)
                    {
                        mask = cValue << (bitsRemaining - 5);
                        curByte = (byte)(curByte | mask);
                        bitsRemaining -= 5;
                    }
                    else
                    {
                        mask = cValue >> (5 - bitsRemaining);
                        curByte = (byte)(curByte | mask);
                        returnArray[arrayIndex++] = curByte;
                        curByte = (byte)(cValue << (3 + bitsRemaining));
                        bitsRemaining += 3;
                    }
                }

                //if we didn't end with a full byte
                if (arrayIndex != byteCount)
                {
                    returnArray[arrayIndex] = curByte;
                }

                return returnArray;
            }

            public static string ToString(byte[] input)
            {
                if (input == null || input.Length == 0)
                {
                    throw new ArgumentNullException("input");
                }

                int charCount = (int)Math.Ceiling(input.Length / 5d) * 8;
                char[] returnArray = new char[charCount];

                byte nextChar = 0, bitsRemaining = 5;
                int arrayIndex = 0;

                foreach (byte b in input)
                {
                    nextChar = (byte)(nextChar | (b >> (8 - bitsRemaining)));
                    returnArray[arrayIndex++] = ValueToChar(nextChar);

                    if (bitsRemaining < 4)
                    {
                        nextChar = (byte)((b >> (3 - bitsRemaining)) & 31);
                        returnArray[arrayIndex++] = ValueToChar(nextChar);
                        bitsRemaining += 5;
                    }

                    bitsRemaining -= 3;
                    nextChar = (byte)((b << bitsRemaining) & 31);
                }

                //if we didn't end with a full char
                if (arrayIndex != charCount)
                {
                    returnArray[arrayIndex++] = ValueToChar(nextChar);
                    while (arrayIndex != charCount) returnArray[arrayIndex++] = '='; //padding
                }

                return new string(returnArray);
            }

            private static int CharToValue(char c)
            {
                int value = (int)c;

                //65-90 == uppercase letters
                if (value < 91 && value > 64)
                {
                    return value - 65;
                }
                //50-55 == numbers 2-7
                if (value < 56 && value > 49)
                {
                    return value - 24;
                }
                //97-122 == lowercase letters
                if (value < 123 && value > 96)
                {
                    return value - 97;
                }

                throw new ArgumentException("Character is not a Base32 character.", "c");
            }

            private static char ValueToChar(byte b)
            {
                if (b < 26)
                {
                    return (char)(b + 65);
                }

                if (b < 32)
                {
                    return (char)(b + 24);
                }

                throw new ArgumentException("Byte is not a value Base32 value.", "b");
            }

        }
        [STAThread]
        static void Main(string[] args)
        {
            var bytes = Base32Encoding.ToBytes("JBSWY3DPEHPK3PXP");

            var totp = new Totp(bytes);

            var result = totp.ComputeTotp();
            Clipboard.SetText(result.ToString());

            //var remainingTime = totp.RemainingSeconds();
        }
    }
}