damyarou

python, GMT などのプログラム

Python jpg画像を縮小保存しhtmlで表示

タイトル通り、jpg画像を縮小保存しhtmlで表示する。

import glob
import os
from PIL import Image
import os

def select_pic():
    path = '*.jpg'
    # 現フォルダ内にあって拡張子がJPGのファイル名を取得
    file_list = glob.glob(path)
    print(file_list)
    return file_list

def process(file_list):
    print(file_list)
    #for frameR in file_list:
    for i in range(len(file_list)):
        fnameR=file_list[i]
        print(fnameR)
        # 元画像の読み込み
        originalImg = Image.open(fnameR)
        # 元画像のサイズを取得
        width1, height1 = originalImg.size
        # サムネイル作成
        width2=600
        height2=width2*height1/width1
        originalImg.thumbnail((int(width2), int(height2)), Image.ANTIALIAS)
        name=fnameR.replace('.jpg','')
        fnameW='600_'+name+'.jpg'
        # 保存
        originalImg.save(fnameW)


def make_html():
    filenames = os.listdir('./')
    imgl=[]
    ww=[]
    hh=[]
    for fname in sorted(filenames):
        if fname[0:4]=='600_':
            print(fname)
            im=Image.open(fname)
            w=im.size[0]
            h=im.size[1]
            imgl=imgl+[fname]
            ww=ww+[w]
            hh=hh+[h]
    f=open('maggie.html','w')
    print('<html>',file=f)
    print('<body>',file=f)
    print('<table>',file=f)
    n=len(imgl)
    m=int(n/5)+1
    k=-1
    for i in range(0,m):
        print('<tr>',file=f)
        for j in range(0,5):
            k=k+1
            if k<=n-1:
                pic=imgl[k]
                w1=200
                h1=int(hh[k]/ww[k]*200)
                print('<td align="center"><img src="'+pic+'" alt="pic" width="'+str(w1)+'", height="'+str(h1)+'"><br><a href="'+pic+'">'+pic+'<a></td>',file=f)
            else:
                print('<td></td>',file=f)
        print('</tr>',file=f)
    print('</table>',file=f)
    print('</body>',file=f)
    print('</html>',file=f)
    f.close()

def main():
    file_list=select_pic()
    process(file_list)
    make_html()


if __name__ == '__main__':
    main()

以　上

雑記　また機械学習の本を買ってしまった

記事の最後に行く

昨日（2019.08.11）、また機械学習の本を買ってしまった。

あと、久しぶりにpipで一括アップデートしたら、matplotlibがなんとなくおかしい。

前にjupyterがおかしかったときは、tornadoの問題だった（以下参照）

今回はmatplotlibの問題のようです。

Thank you.

記事の先頭に行く

雑記 機械学習の勉強を始めてみる

記事の最後に行く

いよいよ７月。 昨年（2018年10月10日）、マレーシアでの業務が完了し、帰国した。 それ以降、以下のような感じで過ごしてきた。

  • 2018年10月 ひま
  • 2018年11月 ひま
  • 2018年12月 忙しい
  • 2019年01月 忙しい
  • 2019年02月 忙しい（ベトナム出張２週間）
  • 2019年03月 忙しい
  • 2019年04月 上旬は忙しい（マレーシア出張２週間）、中旬以降ひま
  • 2019年05月 ひま
  • 2019年06月 ひま

昨年の10月11月は急速の意味もありひまでよかったのだが、今年になって4月中旬以降の2ヶ月半はなんだか時間を無駄に過ごしてきた感じ。

そこで昨日久しぶりに本を買った。

scikit-learn を用いて機械学習の基礎を勉強するもの。

機械学習もちらほらネットを見ながらかじっては見ていたが、やはり本に沿って一式流してみないと体系的な知識がつかないなと感じていたところ。

今の予定では８月からまた忙しくなる「予定」ではあるが、第１章を読んだところ、比較的わかりやすく書かれているようなので、なんとか最後まで読破しようと思っている。

ついでに英語の勉強も少しやってみようかな。

Thank you.

記事の先頭に行く

設計　Pythonによる水文量頻度解析用自作関数とその利用

設計

記事の最後に行く

自作関数

以下の確率分布に対するパラメータ推定と予測値計算の部分を自作関数としてプログラム本体から分離している。

含まれる関数（確率分布）
関数名説明
ln3３変数対数正規分布（LN3）
lp3対数ピアソンIII型分布（LP3）
gev一般化極値分布（GEV）
gumグンベル分布（GUM）
sqe平方根指数型最大値分布（SQRT-ET）

関数への入力値（引数）はいずれについても以下の３つ。

入力値（引数）
変数名説明
xx昇順に並び替えられた水文量観測値
pp各観測値に相当する非超過確率（プロッティング・ポジション公式より算出）<.td>
pyear予測値を求めたい確率年（昇順）

関数からの戻り値はいずれについても以下の４つ。

戻り値
変数名説明
res計算された確率分布関数のパラメータ（numpy.array）
ye入力した非超過確率（pp）に対応する予測値（numpy.array）
yp入力した確率年に対する予測値（numpy.array）
ss確率分布名、算定されたパラメータなどの情報を記載した文字列


#im_hydrology.py
import numpy as np
from scipy.stats import norm  # normal distribution for LN3
from scipy.stats import gamma # Gamma distribution for LP3
import scipy.special as ssp   # Gamma function for GEV

def ln3(xx,pp,pyear):
    def func_ln3(xx):
        n=len(xx)
        mx=np.sum(xx)/n
        s1=np.sqrt(np.sum((xx-mx)**2)/n)
        cx=np.sum(((xx-mx)/s1)**3)/n
        sx=np.sqrt(n/(n-1))*s1
        #rx=np.sqrt(n*(n-1))/(n-2)*cx
        a=1.01+7.01/n+14.66/n**2
        b=1.69/n+74.66/n**2
        rx=cx*(a+b*cx**3)
        beta=1+rx**2/2
        x=(beta+np.sqrt(beta**2-1))**(1/3)+(beta-np.sqrt(beta**2-1))**(1/3)-1
        sd=np.sqrt(np.log(x))
        mu=np.log(sx/np.sqrt(x*(x-1)))
        aa=mx-np.exp(mu)*np.exp(sd**2/2)
        return aa,mu,sd
    def est_ln3(aa,mu,sd,pp):
        zz=norm.ppf(pp, loc=0, scale=1)
        return aa+np.exp(mu+sd*zz)
    # parameter calculation
    aa,mu,sd=func_ln3(xx)
    # estimation
    ye=est_ln3(aa,mu,sd,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_ln3(aa,mu,sd,prob)
    res=np.array([aa,mu,sd,rr])
    ss='LN3\na={0:8.3f}\nm={1:8.3f}\ns={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def lp3(xx,pp,pyear):
    def func_lp3(xx):
        n=len(xx)
        yy=np.log(xx)
        my=np.sum(yy)/n
        sy=np.sqrt(np.sum((yy-my)**2)/n)
        cy=np.sum(((yy-my)/sy)**3)/n
        sy=np.sqrt(n/(n-1))*sy
        a=1+6.51/n+20.2/n**2
        b=1.48/n+6.77/n**2
        ry=cy*(a+b*cy**2)
        bb=4/ry**2
        aa=sy/np.sqrt(bb)
        if ry<0: aa=-aa
        cc=my-aa*bb
        return bb,cc,aa
    def est_lp3(bb,cc,aa,pp):
        ww=gamma.ppf(pp, bb, loc=0, scale=1)
        return np.exp(cc+aa*ww)
    # parameter calculation
    bb,cc,aa=func_lp3(xx)
    # estimation
    ye=est_lp3(bb,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_lp3(bb,cc,aa,prob)
    res=np.array([bb,cc,aa,rr])
    ss='LP3\nb={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gev(xx,pp,pyear):
    def func_gev(xx):
        n=len(xx)
        jj=np.arange(n)+1
        b0=np.sum(xx)/n
        b1=np.sum((jj-1)*xx)/n/(n-1)
        b2=np.sum((jj-1)*(jj-2)*xx)/n/(n-1)/(n-2)
        lam1=b0
        lam2=2*b1-b0
        lam3=6*b2-6*b1+b0
        d=2*lam2/(lam3+3*lam2)-np.log(2)/np.log(3)
        k=7.8590*d+2.9554*d**2
        a=k*lam2/(1-2**(-k))/ssp.gamma(1+k)
        c=lam1-a/k*(1-ssp.gamma(1+k))
        return k,c,a
    def est_gev(kk,cc,aa,pp):
        return cc+aa/kk*(1-(-np.log(pp))**kk)
    # parameter calculation
    kk,cc,aa=func_gev(xx)
    # estimation
    ye=est_gev(kk,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gev(kk,cc,aa,prob)
    res=np.array([kk,cc,aa,rr])
    ss='GEV\nk={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gum(xx,pp,pyear):
    def func_gum(xx):
        n=len(xx)
        b0=np.sum(xx)/n
        j=np.arange(0,n)
        b1=np.sum(j*xx)/n/(n-1)
        lam1=b0
        lam2=2*b1-b0
        aa=lam2/np.log(2)
        cc=lam1-0.5772*aa
        return aa,cc
    def est_gum(aa,cc,pp):
        return cc-aa*np.log(-np.log(pp))
    aa,cc=func_gum(xx)
    # estimation
    ye=est_gum(aa,cc,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gum(aa,cc,prob)
    res=np.array([aa,cc,rr])
    ss='GUM\na={0:8.3f}\nc={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss


def sqe(xx,pp,pyear):
    def est_sqrtet(aa,bb,pp,x0):
        # estimation of quantile by Newton-Raphson Method
        def gx(x,aa,p):
            ggx=np.log(1+x)-x-np.log(-1/aa*np.log(p))
            gdx=1/(1+x)-1
            return ggx,gdx
        imax=100
        eps=1.0e-10
        tp=np.zeros(len(pp),dtype=np.float64)
        x1=x0
        for i,p in enumerate(pp):
            for j in range(imax):
                ggx,gdx=gx(x1,aa,p)
                x2=x1-ggx/gdx
                if np.abs(x2-x1)<eps: break
                x1=x2
            tp[i]=x2
        return tp**2/bb

    def bisection(xx):
        # estimation of coefficient bb by Bisection Method
        def func(bb,xx):
            # function for Bisection Method
            nn=len(xx)
            a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
            a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx))))
            return a1-a2
        imax=100    # maximum number of iterations
        err=1.0e-10 # allowable error of the zero value
        x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2
        x2=x1+0.5
        for i in range(0,imax):
            xi=0.5*(x1+x2)
            f1=func(x1,xx)
            f2=func(x2,xx)
            fi=func(xi,xx)
            if 0<f1*f2:
                x1=x2
                x2=x1+0.5
            elif abs(x2-x1) < err:
                break # converged
            else:
                if f1*fi < 0.0: x2=xi
                if fi*f2 < 0.0: x1=xi
        return xi

    bb=bisection(xx)
    aa=(np.sum(np.sqrt(bb*xx))-2*len(xx))/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
    # estimation
    x0=np.sqrt(bb*xx[-1])
    ye=est_sqrtet(aa,bb,pp,x0)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_sqrtet(aa,bb,prob,x0)
    res=np.array([aa,bb,rr])
    ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss

クオンタイルプロット

f:id:damyarou:20190623235303j:plain

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology


fsz=12
xmin=0; xmax=5000; dx=1000
ymin=0; ymax=5000; dy=1000


def qqplot(fnameR,xx,ye,ss,n):
    tstr=fnameR
    nnn=230+n
    plt.subplot(nnn)
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Observed (m$^3$/s)')
    plt.ylabel('Predicted (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.plot(xx,ye,'o',clip_on=False)
    plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1)
    xs=xmin+0.05*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished')


def calc(fnameR,alpha,pyear):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    res1,ye1,yp1,ss1=im_hydrology.ln3(xx,pp,pyear)    
    res2,ye2,yp2,ss2=im_hydrology.lp3(xx,pp,pyear)    
    res3,ye3,yp3,ss3=im_hydrology.gev(xx,pp,pyear)
    res4,ye4,yp4,ss4=im_hydrology.gum(xx,pp,pyear)
    res5,ye5,yp5,ss5=im_hydrology.sqe(xx,pp,pyear)
    
    plt.figure(figsize=(12,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    qqplot(fnameR,xx,ye1,ss1,1)
    qqplot(fnameR,xx,ye2,ss2,2)
    qqplot(fnameR,xx,ye3,ss3,3)
    qqplot(fnameR,xx,ye4,ss4,4)
    qqplot(fnameR,xx,ye5,ss5,5)
    plt.tight_layout()
    s=fnameR.replace('inp_','')
    s=s.replace('.txt','')
    fnameF='fig_qq_'+s+'.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()    

    s0='{0:8s}'.format('Year')
    s1='{0:8s}'.format('LN3')
    s2='{0:8s}'.format('LP3')
    s3='{0:8s}'.format('GEV')
    s4='{0:8s}'.format('GUM')
    s5='{0:8s}'.format('SQRT-ET')
    for i in range(len(pyear)):
        s0=s0+',{0:6.0f}'.format(pyear[i])
        s1=s1+',{0:6.0f}'.format(yp1[i])
        s2=s2+',{0:6.0f}'.format(yp2[i])
        s3=s3+',{0:6.0f}'.format(yp3[i])
        s4=s4+',{0:6.0f}'.format(yp4[i])
        s5=s5+',{0:6.0f}'.format(yp5[i])
    print(s0)
    print(s1)
    print(s2)
    print(s3)
    print(s4)
    print(s5)


def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        calc(fnameR,alpha,pyear)


#==============
# Execution
#==============
if __name__ == '__main__': main()

確率年と予測値

f:id:damyarou:20190623235336j:plain

# Jackknife method
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology

# Global variables
fsz=14
xmin=1
xmax=1000
dx=1
ymin=0
ymax=10000
dy=1000


def drawfig0(pyear,qqq,smt,ss):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qqq[0,:],'o--',ms=10,color=col,markerfacecolor=col,label=smt[0],lw=1,clip_on=False)
    plt.plot(pyear,qqq[1,:],'s--',ms=10,color=col,markerfacecolor=col,label=smt[1],lw=1,clip_on=False)
    plt.plot(pyear,qqq[2,:],'^--',ms=10,color=col,markerfacecolor=col,label=smt[2],lw=1,clip_on=False)
    plt.plot(pyear,qqq[3,:],'o-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[3],lw=1,clip_on=False)
    plt.plot(pyear,qqq[4,:],'s-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[4],lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    fnameF='fig_comp.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def drawfig(pyear,qb,sd,ss,kkk):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    if kkk==1: tstr=tstr+' (LN3)'
    if kkk==2: tstr=tstr+' (LP3)'
    if kkk==3: tstr=tstr+' (GEV)'
    if kkk==4: tstr=tstr+' (GUM)'
    if kkk==5: tstr=tstr+' (SQRT-ET)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qb   ,'o-'  ,ms=10,color=col,markerfacecolor=col,label='Average',lw=1,clip_on=False)
    plt.plot(pyear,qb+sd,'^--' ,ms=10,color=col,markerfacecolor=col,label='Ave.+ SD',lw=1,clip_on=False)
    plt.plot(pyear,qb-sd,'s--' ,ms=10,color=col,markerfacecolor=col,label='Ave.- SD',lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    if kkk==1: fnameF='fig_rp_ln3.jpg'
    if kkk==2: fnameF='fig_rp_lp3.jpg'
    if kkk==3: fnameF='fig_rp_gev.jpg'
    if kkk==4: fnameF='fig_rp_gum.jpg'
    if kkk==5: fnameF='fig_rp_sqe.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def calc(fnameR,alpha,pyear,kkk):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    if kkk==1: res,ye0,yp0,ss=im_hydrology.ln3(xx,pp,pyear)    
    if kkk==2: res,ye0,yp0,ss=im_hydrology.lp3(xx,pp,pyear)    
    if kkk==3: res,ye0,yp0,ss=im_hydrology.gev(xx,pp,pyear)    
    if kkk==4: res,ye0,yp0,ss=im_hydrology.gum(xx,pp,pyear)
    if kkk==5: res,ye0,yp0,ss=im_hydrology.sqe(xx,pp,pyear)
    print(ss)
    
    qq=np.zeros((n,len(pyear)),dtype=np.float64)
    qm=np.zeros(len(pyear),dtype=np.float64)
    qb=np.zeros(len(pyear),dtype=np.float64)
    sd=np.zeros(len(pyear),dtype=np.float64)
    for i in range(0,n):
        xi=np.delete(xx,i)
        m=len(xi)
        ii=np.arange(m)+1
        pp=(ii-alpha)/(m+1-2*alpha)
        if kkk==1: res,ye,yp,ss=im_hydrology.ln3(xi,pp,pyear)
        if kkk==2: res,ye,yp,ss=im_hydrology.lp3(xi,pp,pyear)
        if kkk==3: res,ye,yp,ss=im_hydrology.gev(xi,pp,pyear)
        if kkk==4: res,ye,yp,ss=im_hydrology.gum(xi,pp,pyear)
        if kkk==5: res,ye,yp,ss=im_hydrology.sqe(xi,pp,pyear)
        qq[i,:]=yp
    for j in range(0,len(pyear)):
        qm[j]=np.average(qq[:,j])
        qb[j]=n*yp0[j]-(n-1)*qm[j]
        sd[j]=np.sqrt((n-1)/n*np.sum((qq[:,j]-qm[j])**2))
    return qb,sd
    

def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    qqq=np.zeros((5,len(pyear)),dtype=np.float64)
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        for kkk in [1,2,3,4,5]:
            qb,sd=calc(fnameR,alpha,pyear,kkk)
            s0='{0:4s}'.format('Year')
            s1='{0:4s}'.format('Ave.')
            s2='{0:4s}'.format('SD')
            for j in range(len(pyear)):
                s0=s0+'{0:6.0f}'.format(pyear[j])
                s1=s1+'{0:6.0f}'.format(qb[j])
                s2=s2+'{0:6.0f}'.format(sd[j])
            ss='Retuen Period and Flood Discharge (m$^3$/s)'
            ss=ss+'\n'+s0+'\n'+s1+'\n'+s2
            ss=ss+'\n* Ave.: average, SD: standard deviation'
            print(ss)
            drawfig(pyear,qb,sd,ss,kkk)
            qqq[kkk-1,:]=qb[:]
        smt=['LN3','LP3','GEV','GUM','SQRT-ET']
        s0='{0:8s}'.format('Year')
        s1='{0:8s}'.format(smt[0])
        s2='{0:8s}'.format(smt[1])
        s3='{0:8s}'.format(smt[2])
        s4='{0:8s}'.format(smt[3])
        s5='{0:8s}'.format(smt[4])
        for j in range(len(pyear)):
            s0=s0+'{0:6.0f}'.format(pyear[j])
            s1=s1+'{0:6.0f}'.format(qqq[0,j])
            s2=s2+'{0:6.0f}'.format(qqq[1,j])
            s3=s3+'{0:6.0f}'.format(qqq[2,j])
            s4=s4+'{0:6.0f}'.format(qqq[3,j])
            s5=s5+'{0:6.0f}'.format(qqq[4,j])
        ss='Retuen Period and Flood Discharge (m$^3$/s)'
        ss=ss+'\n'+s0+'\n'+s1+'\n'+s2+'\n'+s3+'\n'+s4+'\n'+s5
        print(ss)
        drawfig0(pyear,qqq,smt,ss)


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く

FEM 格子桁構造解析プログラム

FEM

記事の最後に行く

はじめに

格子桁構造解析プログラムを紹介します。 このプログラムは、平面骨組構造解析プログラムの軸力・軸方向変位をねじりモーメント・ねじり角に書き換えたものです。 構造は平面形状とし、x - y 平面上で定義します。

格子桁解析は、立体骨組構造解析プログラムで行うことができるのですが、床あるいは壁のような単純な板状構造物の解析には、入力データの単純さを考えると格子桁解析も有用であることから、作ったものです。

なお、このプログラムでは温度変化による断面力は考慮できません。 それは、考慮できる断面力は、ねじりモーメント、曲げモーメント、せん断力であり、軸力項が無いためです。

（注意事項）

充実矩形断面のねじり定数 J は、断面の短辺を a、長辺を b として、以下の式で計算できます。

\begin{equation} J=\cfrac{1}{3} b a^3 \left(1-\cfrac{2}{\pi} \cfrac{a}{b}\right) \end{equation}

ここで、コンクリート床のような版構造を格子桁でモデル化する場合、ねじり定数 J はゼロとみなせるくらいに小さな値を入力します。 これは、ある要素のねじり角は、それに交差する要素の曲げ角に影響するため、曲げモーメントがねじりの影響で不連続になってしまうためです。

本来の格子桁をモデル化する場合は、もちろん、計算により求めた正規のねじり定数を入力する必要があります。

Grid Girder Analysis

Outline of Grid Girder Analysis Program

  • This is a program for Grid Girder Analysis. Elastic and small displacement problems can be treated. This program has been made using 'Python 3.'
  • Beam element with 2 nodes is used. 1 node has 3 degrees of freedom which are torsional rotation, bending rotation and deflection.
  • Simultaneous linear equations are solved using numpy.linalg.solve.
  • Input / Output file name can be defined arbitrarily as text files with the separator of blank, and those are inputted from command line of a terminal.
Workable condition
ItemDescription
ElementBeam element which has 3 degrees of freedom such as torsional rotation, bending rotation and deflection.
LoadSpecify the node number and load values on the node.
DisplacementSpecify the node number and displacements on the node. Any values can be applied including zero


f:id:damyarou:20190605111530p:plain f:id:damyarou:20190605111546p:plain
Global Coordinate SystemLocal Coordinate Syatem


The Stiffness matrix of an element is shown below:

\begin{equation} \begin{Bmatrix} T_i \\ M_i \\ Q_i \\ T_j \\ M_j \\ Q_j \end{Bmatrix} =\begin{bmatrix} GJ/L & 0 & 0 & -GJ/L & 0 & 0 \\ 0 & 4 EI/L & -6 EI/L^2 & 0 & 2 EI/L & 6 EI/L^2 \\ 0 & -6 EI/L^2 & 12 EI/L^3 & 0 & -6 EI/L^2 & -12 EI/L^3 \\ -GJ/L & 0 & 0 & GJ/L & 0 & 0 \\ 0 & 2 EI/L & -6 EI/L^2 & 0 & 4 EI/L & 6 EI/L^2 \\ 0 & 6 EI/L^2 & -12 EI/L^3 & 0 & 6 EI/L^2 & 12 EI/L^3 \end{bmatrix} \begin{Bmatrix} \phi_i \\ \theta_i \\ w_i \\ \phi_j \\ \theta_j \\ w_j \end{Bmatrix} \end{equation}


\begin{align} &GJ & &\text{: Torsional rigidity} & &T & &\text{: Torsional moment} & &\phi & &\text{: Tortional rotation around x-axis}\\ &EI & &\text{: Bending rigidity} & &M & &\text{: Bending moment} & &\theta & &\text{: Bending rotation around y-axis}\\ &L & &\text{: Length of element} & &Q && \text{: Shearing force} & &w & &\text{: Displacement in z-direction} \end{align}

A shear modulus G is calculated from an elastic modulus of E and Poisson's ration \nu using following equation in the program.

\begin{equation} G=\cfrac{E}{2 (1+\nu)} \end{equation}

Since the coordinate transformation is carried out on the x-y plane only, the transfoemation matrix is the same type as it for 2D frame analysis.

Format of input data file

npoin nele nsec npfix nlod  # Basic values for analysis
E po AA AI AJ gamma         # material properties
    ....(1 to nsec)....     #
node-1,node-2,matno         # Element connectivity, material set number, uniformly distributed load per unit length
    ....(1 to nele)....     #
x,y                         # Node coordinates, temperature change of node
    ....(1 to npoin)....    #
node nokx noky nokz rdisx rdisy rdisz # node number, restrict conditions and forced displacement
    ....(1 to npfix)....    #
node,fx,fy,fz               # node number, Load values of torsional moment, bending moment, vertical force
    ....(1 to nlod)....     # (Omit data input if nlod=0)


npoin : Number of nodes                   E     : Elastic modulus of element
nele  : Number of elements                po    : Poison's ratio
nsec  : Number of material sets           AA    : Section area of element
npfix : Number of restricted nodes        AI    : Moment of Inertia of element
nlod  : Number of loaded nodes            AJ    : Tosion constant
                                          gamma : unit weight of material
                                          matno : Material set number
  • The structure shall be defined on the x-y plane in the global coordinate system.
  • x-direction in the local coordinate system means the member axis.
  • The displacement in local x-direction means the rotation around the x-axis, which is same as torsional rotation.
  • The displacement in local y-direction means the rotation around the y-axis, which is same as bending rotation.
  • The displacement in local z-direction means the deflection of the beam.
  • Restricted node means the node which has known (given) displacement. As a known (given) value of nodal displacement, any value can be given including zero for a restricted node.
  • Since stress resultants of element are defined as equivalent nodal forces in local coordinate system, it is necessary to note that signs are different from it on general structural mechanics. Positive directions are always right-direction, upward-direction for each node in local coordinate system.

Format of output file

npoin  nele  nsec npfix  nlod
    (Each value for above item)
  sec   E   po   A   I   J   gamma
    (material caracteristics)
    .....(1 to nsec).....
 node   x   y   fx   fy   fz   kox   koy   koz
    node: node number
    x   : x-coordinate
    y   : y-coordinate
    fx  : torsional load
    fy  : bending moment load
    fz  : vertical force load
    kox : restriction around x-axis
    koy : restriction around y-axis
    koz : restriction in z-direction
    .....(1 to npoin).....
 node   kox   koy   koz   rdis_x   rdis_y   rdis_z
    node   : node number
    kox    : restriction around x-axis
    koy    : restriction around y-axis
    koz    : restriction in z-direction
    rdis_x : given rotation around x-axis
    rdis_y : given rotation around y-axis
    rdis_z : given displacement in z-direction
    .....(1 to npfix).....
 elem   i   j   sec
    elem : element number
    i    : node-1
    j    : node-2
    sec  : material number
    .....(1 to nele).....
 node   rot-x   rot-y   dis-z
    node  : node number
    rot-x : rotation arouns x-axis
    rot-y : rotation around y-axis
    dis-z : displacement in z-direction
    .....(1 to npoin).....
 elem   T_i   M_i   Q_i   T_j   M_j   Q_j
    elem : element number
    T_i  : torsional moment at node-i
    M_i  : bending moment at node-i
    Q_i  : shearing force at node-i
    T_j  : torsional moment at node-j
    M_j  : bending moment at node-j
    Q_j  : shearing force at node-j
    .....(1 to nele).....
n=693  time=0.078 sec   # reference information

Grid Girder Analysis Program

###################################
# Grid Girder Analysis
###################################
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_grd(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    # array declaration
    x     =np.zeros([2,npoin],dtype=np.float64)         # Coordinates of nodes
    ae    =np.zeros([6,nsec],dtype=np.float64)          # Section characteristics
    node  =np.zeros([nod+1,nele],dtype=np.int)          # Node-element relationship
    fp    =np.zeros(nfree*npoin,dtype=np.float64)       # External force vector
    mpfix =np.zeros([nfree,npoin],dtype=np.int)         # Ristrict conditions
    rdis  =np.zeros([nfree,npoin],dtype=np.float64)     # Ristricted displacement
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) # E  : Elastic modulus
        ae[1,i]=float(text[1]) # po : Poisson's ratio
        ae[2,i]=float(text[2]) # AA : Section area
        ae[3,i]=float(text[3]) # AI : Moment of inertia
        ae[4,i]=float(text[4]) # AJ : Tortional constant
        ae[5,i]=float(text[5]) # gamma : Unit weight of beam material
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]=float(text[0])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])             #fixed node
        mpfix[0,lp-1]=int(text[1])  #index of restriction for tortional rotation
        mpfix[1,lp-1]=int(text[2])  #index of restriction for bending rotation
        mpfix[2,lp-1]=int(text[3])  #index of restriction for vertical displacement
        rdis[0,lp-1]=float(text[4]) #given tortional rotation
        rdis[1,lp-1]=float(text[5]) #given bending rotation
        rdis[2,lp-1]=float(text[6]) #given vertical displacement
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[3*lp-3]=float(text[1]) #tortional moment load
            fp[3*lp-2]=float(text[2]) #bending moment load
            fp[3*lp-1]=float(text[3]) #vertical load
    f.close()
    return npoin, nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp


def prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('sec','E','po','A','I','J','gamma'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
    .format('node','x','y','fx','fy','fz','kox','koy','koz'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d} {8:5d}'
        .format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>15s} {5:>15s} {6:>15s}'
    .format('node','kox','koy','koz','rdis_x','rdis_y','rdis_z'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]+mpfix[2,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:15.7e} {5:15.7e} {6:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],rdis[0,i],rdis[1,i],rdis[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

    
def prout_grd(fnameW,npoin,nele,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','rot-x','rot-y','dis-z'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','T_i','M_i','Q_i','T_j','M_j','Q_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
    fout.close()


def sm_grd(GJ,EI,x1,y1,x2,y2):
    ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
    tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
    xx=x2-x1
    yy=y2-y1
    el=np.sqrt(xx**2+yy**2)
    ek[0,0]= GJ/el; ek[0,1]=        0.0; ek[0,2]=         0.0; ek[0,3]=-GJ/el; ek[0,4]=        0.0; ek[0,5]=         0.0
    ek[1,0]=   0.0; ek[1,1]= 4*EI/el   ; ek[1,2]= -6*EI/el**2; ek[1,3]=   0.0; ek[1,4]= 2*EI/el   ; ek[1,5]=  6*EI/el**2
    ek[2,0]=   0.0; ek[2,1]=-6*EI/el**2; ek[2,2]= 12*EI/el**3; ek[2,3]=   0.0; ek[2,4]=-6*EI/el**2; ek[2,5]=-12*EI/el**3
    ek[3,0]=-GJ/el; ek[3,1]=        0.0; ek[3,2]=         0.0; ek[3,3]= GJ/el; ek[3,4]=        0.0; ek[3,5]=         0.0
    ek[4,0]=   0.0; ek[4,1]= 2*EI/el   ; ek[4,2]= -6*EI/el**2; ek[4,3]=   0.0; ek[4,4]= 4*EI/el   ; ek[4,5]=  6*EI/el**2
    ek[5,0]=   0.0; ek[5,1]= 6*EI/el**2; ek[5,2]=-12*EI/el**3; ek[5,3]=   0.0; ek[5,4]= 6*EI/el**2; ek[5,5]= 12*EI/el**3
    s=yy/el
    c=xx/el
    tt[0,0]=  c; tt[0,1]=  s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
    tt[1,0]= -s; tt[1,1]=  c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
    tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
    tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]=  c; tt[3,4]=  s; tt[3,5]=0.0
    tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]=  c; tt[4,5]=0.0
    tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
    return ek,tt


def bfvec_grd(gamma,A,x1,y1,x2,y2):
    # Distributed vertical load vector
    bfv=np.zeros(6,dtype=np.float64)
    el=np.sqrt((x2-x1)**2+(y2-y1)**2)
    gkv=-1.0 # downward direction
    bfv[0]=0.0
    bfv[1]=0.0
    bfv[2]=0.5*gamma*A*el*gkv
    bfv[3]=0.0
    bfv[4]=0.0
    bfv[5]=0.5*gamma*A*el*gkv
    return bfv


def main():
    # Main routine
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=3            # Degree of freedom per node
    # data input and input data print
    npoin,nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp=inpdata_grd(fnameR,nod,nfree)
    prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node)
    # assembly of stiffness matrix & load vectors
    ir    =np.zeros(nod*nfree,dtype=np.int)             # Work vector for matrix assembly
    gk    =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64) # Global stiffness matrix
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        aa=ae[2,m]
        ai=ae[3,m]
        aj=ae[4,m]
        gamma=ae[5,m]
        A=aa
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)      # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)     # Stiffness matrix in global coordinate
        bfe  =bfvec_grd(gamma,A,x1,y1,x2,y2) # Body force vector downward direction
        ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
        ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+ck[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of section force
    fsec=np.zeros([nod*nfree,nele],dtype=np.float64) # Section force vector
    dis =np.zeros(nod*nfree,dtype=np.float64) # work vector for section force calculation
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        ai=ae[3,m]
        aj=ae[4,m]
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)
        dis[0]=disg[3*i]; dis[1]=disg[3*i+1]; dis[2]=disg[3*i+2]
        dis[3]=disg[3*j]; dis[4]=disg[3*j+1]; dis[5]=disg[3*j+2]
        fsec[:,ne]=np.dot(ek,np.dot(tt,dis))
    # print out of result
    prout_grd(fnameW,npoin,nele,disg,fsec)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く

設計 RC長方形断面の曲げ耐荷力計算

設計

記事の最後に行く

はじめに

マレーシア時代に使っていた、RC長方形断面の曲げ耐荷力計算プログラムを掲載しています。 部分安全係数、コンクリートおよび鉄筋の物性は　BS 8110-1 に準拠したものとしています。

日本のコンクリート標準示方書での計算も、部分安全係数と物性値を変えればできるはずなので、後ほど試してみようと思います。

Ultimate Limit Strength Calculation of Rectangular RC member

Stress-strain curve

Stress-strain curves of concrete and rebar follow BS 8110-1.

The compressive stress and strain have positive sign for both of concrete and rebar.

Concrete

\begin{equation} \sigma_c= \begin{cases} \quad 0 & (\epsilon <0) \\ \quad \sigma_{cy}-\cfrac{\sigma_{cy}}{\epsilon_{cy}^2}(\epsilon - \epsilon_{cy})^2 & (0\leqq \epsilon \leqq \epsilon_{cy}) \\ \quad \sigma_{cy} & (\epsilon_{cy} < \epsilon \leqq \epsilon_{cu}) \\ \quad 0 & (\epsilon_{cu} < \epsilon) \end{cases} \end{equation}
\begin{gather} \sigma_{cy}=0.67\cdot \cfrac{f_{cu}}{\gamma_m} \qquad \epsilon_{cy}=2.4\times 10^{-4} \sqrt{\cfrac{f_{cu}}{\gamma_m}} \qquad \gamma_m=1.5 \qquad \epsilon_{cu}=0.0035 \end{gather}

where, f_{cu} means a compressive strength of concrete based on cubic specimen following BS. For example, C35 in BS means f_{cu}=35 (N/mm ^2).

Reinforcement

\begin{equation} \sigma_s= \begin{cases} \quad -\sigma_{sy} & (\epsilon < -\epsilon_{sy}) \\ \quad E_s \cdot \epsilon & (-\epsilon_{sy} \leqq \epsilon \leqq \epsilon_{sy}) \\ \quad \sigma_{sy} & (\epsilon_{sy} < \epsilon) \\ \end{cases} \end{equation}
\begin{gather} \sigma_{sy}=\cfrac{f_y}{\gamma_m} \qquad \epsilon_{sy}=\cfrac{f_y}{\gamma_m E_s} \qquad \gamma_m=1.05 \qquad E_s=2.0\times 10^5 \text{(N/mm$^2$)} \\ \end{gather}

where, f_y means a yield strength of rebar and usually f_y=460 (N/mm ^2) is used.

Stress-Strain Curves of Materials
f:id:damyarou:20190602170951j:plain

Calculation of Axial Force and Bending Moment

Calculation formulas for axial force and bending moment are shown below. The moment is calculated with center point of the section.

\begin{align} &\text{Axial force} &\qquad &N=\int_{0}^{h} \sigma \cdot dA=\sum_{i=1}^{m} \sigma_i \cdot b \cdot dy_i \\ &\text{Bending moment} &\qquad &M=\int_{0}^{h} \sigma \cdot (y-h/2) \cdot dA=\sum_{i=1}^{m} \sigma_i \cdot b \cdot dy_i \cdot (y_i-h/2) \end{align}


Model of Cross Section
f:id:damyarou:20190602171009j:plain:w300

Simplified Ultimate Limit State assuming a strain distribution in a section

At the ultimate limit state of a section, either compressive strain of concrete or tensile strain of rebar reachs its limit ultimate strain.

Considering above, simplified strain distribution at ultimate limit state of a section is assumed shown as above figure.

The calculation process is shown below:

  • Set ultimate limit compressive strain of concrete as \epsilon_{cu}
  • Set ultimate limit tensile strain of rebar as -3 \epsilon_{sy} (conveniently assumed value)
  • Step-1: Fix the upper edge strain \epsilon_1=\epsilon_{cu}. Vary the lower edge strain \epsilon_2 from \epsilon_{cu} to -3 \epsilon_{sy} and calculate the axial force and bending moment.
  • Step-2: Fix the lower edge strain \epsilon_2=-3 \epsilon_{sy}. Vary the upper edge strain \epsilon_1 from \epsilon_{cu} to -3 \epsilon_{sy} and calculate the axial force and bending moment.

Obtained axial forces and bending moments from above process deam the N-M relationship at the ultimate limit state of a section.


Concept of Calculation of Ultimate Strength
f:id:damyarou:20190602171711j:plain

Test Running

Used properties of RC cross section and parameters are shown below.

b=1000 (mm) section width
h=1000 (mm) section height
dt=100 (mm) cover of tensile rebar (from tensile edge to rebar center)
dc=100 (mm) cover of compressive rebar (from compressive edge to rebar center)
p1 (%)ratio of tensile reinforcement sres to RC cross section area
p2 (%)ratio of compressive reinforcement sres to RC cross section area
Concrete gradeC25
Yield strength of reinforcement460 MPa

In the following program, the strain difference between $\epsilon{cu}$ and $-3 \epsilon{sy}$ is divided by 20 (nn=20).

f:id:damyarou:20190602175350j:plain

プログラム

自作関数のインポート

ここでは、曲げ耐荷力を計算する自作関数「im_rebar.py」を作成し、メインプログラムから呼ぶようにしています。 Jupyter 上で作業する場合、計算ケースが変わる場合は違うセルを用いますが、同じ部分を何度もコピーして使っていると、意に反して変更してしまい、トラブルの原因にもなるので、この方法を用いています。

自作関数が、Jupyter を実行するのと同じディレクトリにあれば、以下のように関数をインポートすれば使えます。

import im_rebar

自作関数が、Jupyter を実行するのと違うディレクトリにあれば、以下のように関数が格納されているディレクトリのフルパスを import sys で加えた後に関数をインポートします。

import sys; sys.path.append('/Users/kk/pypro')
import im_rebar

負の曲げモーメント

引張鉄筋量と圧縮鉄筋量が等しい場合は負のモーメントは発生しませんが、これらが異なる場合、圧縮応力あるいは引張応力が大きい場合、鉄筋量の差により負のモーメントが発生します。この問題の処理として、このプログラムでは負のモーメントは無視しています。すなわち断面耐荷力の有効範囲は、計算上の正のモーメントの領域のみとしています。

自作関数（im_rebar.py）

# im_rebar.py
import numpy as np


def calmain(mh,fcu,bb,hh,asr,ys):
    # Properties of concrete
    rmc=1.5    # partial safety factor for concrete
    ecu=0.0035 # ultimate limit strain of concrete
    # properties of rebar
    fsy=460.0  # yield strength of rebar
    rms=1.05   # partial safety factor for rebar
    Es=2.0e5   # elastic modulus of steel
    # Division of section height
    yg=hh/2     # location of gravity center
    y=np.zeros(mh,dtype=np.float64)
    dy=hh/float(mh)
    y=np.linspace(0.5*dy,hh-0.5*dy,mh) # distance from tensile edge to each strip
    # Definition of strain for calculation
    ec_lim=1*ecu         # assumed limit compressive strain
    et_lim=-3*fsy/rms/Es # assumed limit tensile strain
    nn=20            # number of division between ec_lim and et_lim
    _e=np.linspace(et_lim,ec_lim,nn+1)
    epsilon=_e[::-1]    # strain values for calculation
    fn=np.zeros(2*(nn+1),dtype=np.float64)  # memory of axial force
    fm=np.zeros(2*(nn+1),dtype=np.float64)  # memory of bending moment
    for (k,ee) in enumerate(epsilon):
        _fn1,_fm1=cal_nm(ec_lim,ee,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es) # Step1 (compressive strain fixed)
        _fn2,_fm2=cal_nm(ee,et_lim,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es) # Step2 (tensile strain fixed)
        fn[k]=_fn1
        fm[k]=_fm1
        fn[k+nn+1]=_fn2
        fm[k+nn+1]=_fm2
    # Arrangement for output
    fnn=fn/bb/hh     # axial force for plotting in stress unit
    fmm=fm/bb/hh**2  # bending moment for plotting in stress unit
    fnc,fnt,fm0=src(nn,fnn,fmm) # search the limit axial forces and pure bending moment
    _i=np.where(0<=fmm)
    fnnp=np.zeros(len(fnn[_i])+2,dtype=np.float64)
    fmmp=np.zeros(len(fmm[_i])+2,dtype=np.float64)
    fnnp[0]=fnc
    fmmp[0]=0
    fnnp[1:len(fnn[_i])+1]=fnn[_i]
    fmmp[1:len(fmm[_i])+1]=fmm[_i]
    fnnp[len(fnn[_i])+1]=fnt
    fmmp[len(fmm[_i])+1]=0
    return fnnp,fmmp,fm0

def sigc(ee,ecu,fcu,rmc):
    # Stress-strain curve of concrete
    scy=0.67*fcu/rmc
    ecy=2.4e-4*np.sqrt(fcu/rmc)
    sig=scy-scy/ecy**2*(ee-ecy)**2
    if ee<=0.0: sig=0.0
    if ecu<=ee: sig=0.0
    if ecy<ee and ee<=ecu: sig=scy
    return sig

def sigs(ee,fsy,rms,Es):
    # Stress-strain curve of Steel rebar
    ssy=fsy/rms
    esy=fsy/rms/Es
    sig=Es*ee
    if ee<-esy: sig=-ssy
    if esy<ee: sig=ssy
    return sig

def cal_nm(e1,e2,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es):
    # Calculation of axial force and bending moment
    fn=0.0
    fm=0.0
    sc=np.zeros(len(y),dtype=np.float64) # concrete stress
    eps=e1-(e1-e2)/hh*(hh-y)
    for i,ee in enumerate(eps):
        sc[i]=sigc(ee,ecu,fcu,rmc)
    ss=np.zeros(len(ys),dtype=np.float64) # rebar stress
    eps=e1-(e1-e2)/hh*(hh-ys)
    for i,ee in enumerate(eps):
        ss[i]=sigs(ee,fsy,rms,Es)-sigc(ee,ecu,fcu,rmc)
    fn=np.sum(sc*dy*bb)+np.sum(ss*asr)
    fm=np.sum(sc*dy*bb*(y-yg))+np.sum(ss*asr*(ys-yg))
    return fn,fm

def src(nn,fnn,fmm):
    # search limit axial forces and pure bending moment by interpolation
    fnc=fnn[0]
    fnt=fnn[-1]
    for k in range(0,2*(nn+1)-1):
        if fmm[k]<0.0 and 0.0<fmm[k+1]:
            fnc=fnn[k]+(fnn[k+1]-fnn[k])/(fmm[k+1]-fmm[k])*(-fmm[k])
        if fmm[k+1]<0.0 and 0.0<fmm[k]:
            fnt=fnn[k]-(fnn[k]-fnn[k+1])/(fmm[k]-fmm[k+1])*(fmm[k])
        if 0.0<fnn[k] and fnn[k+1]<0.0: # estimate ultimate pure bending moment
            fm0=fmm[k]-(fmm[k]-fmm[k+1])/(fnn[k]-fnn[k+1])*fnn[k]
    return fnc,fnt,fm0

メインプログラム

#=============================================
# Ultimate strehgth of Rectangular RC section
# (use module: im_rebar)
#=============================================
import numpy as np
import matplotlib.pyplot as plt
import time
#import sys; sys.path.append('/Users/kk/pypro')
import im_rebar


def sarrow(x1,y1,x2,y2,ss,fsz):
    sv=0
    plt.annotate(
        ss,
        xy=(x1,y1),
        xycoords='data',
        xytext=(x2,y2),
        textcoords='data',
        fontsize=fsz,
        ha='left', va='center',
        arrowprops=dict(shrink=sv,width=0.1,headwidth=5,headlength=10,
                        connectionstyle='arc3',facecolor='#777777',edgecolor='#777777'))
    

def main():
    start=time.time()
    bb=1000
    hh=1000
    mh=int(hh/20)
    fcu=25.0   # design strength of concrete
    lfile=['fig_nm0.jpg','fig_nm1.jpg','fig_nm2.jpg','fig_nm3.jpg']
    lco=['#cccccc','#bbbbbb','#aaaaaa','#999999','#888888','#777777','#666666','#444444','#222222','#000000"']
    pp1=np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]) # tensile bar
    pp2=np.array([0.0, 1.0, 2.0, 3.0]) # compressive bar
    
    fsz=16
    xmin=0; xmax=20; dx=5
    ymin=-40; ymax=40; dy=10
    
    for fnameF,p2 in zip(lfile,pp2):
        tstr='B x H = 1000 x 1000 (C{0:.0f}, p2 = {1:.1f} %)'.format(fcu,p2)
        fig=plt.figure(figsize=(6,6),facecolor='w')
        plt.rcParams["font.size"] = fsz
        plt.rcParams['font.family'] ='sans-serif'
        plt.xlim([xmin,xmax])
        plt.ylim([ymin,ymax])
        plt.xlabel('Bending moment $M\: / \:bh^2$ (MPa)')
        plt.ylabel('Axial force $N\: / \:bh$ (MPa)')
        plt.xticks(np.arange(xmin,xmax+dx,dx))
        plt.yticks(np.arange(ymin,ymax+dy,dy))
        plt.grid(which='major',color='#777777',linestyle='dashed')

        for p1,col in zip(pp1,lco):
            asr=np.array([p1*bb*hh*0.01,p2*bb*hh*0.01])
            ys =np.array([100,hh-100])
            fnp,fmp,fm0=im_rebar.calmain(mh,fcu,bb,hh,asr,ys)
            plt.plot(fmp,fnp,'-',color=col,label='{0:.1f}'.format(p1))
            if p1==0.5:
                plt.fill_betweenx(fnp, 0, fmp, color='#ffffcc')
                ss='Available area\nfor p1=0.5%'
                sarrow(np.mean(fmp),np.mean(fnp),xmin+0.40*(xmax-xmin),ymin+0.85*(ymax-ymin),ss,fsz-2)
            
        plt.legend(loc='lower right',fontsize=fsz-2,title='p1 (%)').get_title().set_fontsize(fsz-2)
        plt.title(tstr,loc='left',fontsize=fsz-2)
        plt.tight_layout()
        plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
        plt.show()    
    print(time.time()-start)


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く

設計 確率雨量の算定（４）追補

設計

記事の最後に行く

平方根指数型最大値分布による確率雨量計算

最近ある方から、昔作った文書http://civilyarou.web.fc2.com/WANtaroHP_html5_win/f90_ENGI/dir_HFA/suimon.pdfの平方根指数型最大値分布に関連して、メールを頂いた。 そのころは、C あるいは Fortran90 を使ってプログラムを作っていたのだが、最近使っている Python で書いてみたくなった。 このため、計算理論は焼き直しであるが、あたらしく Python でプログラムを作ってみた。 ２変数確率分布であるが、未知パラメータを算定するには工夫が必要である。

プログラムとしては、観測データからプロッティングポジションを通さずに直接未知パラメータを計算する従来法と、プロッティングポジションを通して確率分布関数に含まれる未知パラメータを、Scipy の非線形最小二乗法で求める方法の、２種類を紹介している。

計算理論

SQRT exponential-type distribution of maximum (SQRT-ET)

Probabulity Density Function

\begin{equation} f(x)=\cfrac{a b}{2}\cdot \exp\left\{-\sqrt{b x}-a\left(1+\sqrt{b x}\right)\cdot\exp\left(-\sqrt{b x}\right)\right\} \qquad (x \geqq 0) \end{equation}

Cumulative Distribution Function

\begin{equation} F(x)=\exp\left\{-a\left(1+\sqrt{b x}\right)\cdot\exp\left(-\sqrt{b x}\right)\right\} \qquad (x \geqq 0) \end{equation}

Quantile x_p for non-exceedance probability p

\begin{gather} p=\exp\left\{-a\left(1+\sqrt{b x}\right)\cdot\exp\left(-\sqrt{b x}\right)\right\}=\exp\left[-a(1+t)\cdot\exp(-t)\right] \qquad \left(t=\sqrt{b x}\right) \\ \rightarrow \quad x=\cfrac{t{}^2}{b} \qquad \ln(1+t)-t=\ln\left[-\cfrac{1}{a}\ln(p)\right] \end{gather}


\begin{equation} x_p=\cfrac{t_p{}^2}{b} \qquad \ln(1+t_p)-t_p=\ln\left[-\cfrac{1}{a}\ln(p)\right] \end{equation}

In above, the conditions of a \gt 0 and b \gt 0 shall be satisfied, because the terms in logarithmic function and square-root have to be positive values.

Calculation Method of t_p

Since the relationship between t and p is non-linear, Newton-Raphson Method is used to obtain the value of t at non-exceedance probability of p. Following functions are considered, where g' is the differential of function g.

\begin{align} &g(t)=\ln(1+t)-t-\ln\left[-\cfrac{1}{a}\ln(p)\right] \\ &g'(t)=\cfrac{1}{1+t}-1 \end{align}

To obtain the value of t at g(t)=0, iterative calculation will be carried out using following equation.

\begin{equation} t_{i+1}=t_{i}-\cfrac{g(t_{i})}{g'(t_{i})} \qquad i \text{ : counter for iterative calculation} \end{equation}

As an initial value of t, t_0=\sqrt{b\cdot x_{max}} can be used.

Estimation of Parameters (Maximum Likelihood Method)

The parameters a and b are determined with a condition that a log-likelihoog function L shown below has the maximum value.

\begin{align} L(a,b)=&\sum_{j=0}^N \ln f(x_j) \\ =&N\ln a + N\ln b - N\ln 2 - \sum_{j=1}^N \sqrt{b x_j} - a \left\{\sum_{j=1}^N \exp\left(-\sqrt{b x_j}\right) + \sum_{j=1}^N \left[\sqrt{b x_j}\cdot \exp\left(-\sqrt{b x_j}\right)\right]\right\} \end{align}

From the conditions that a partial differential of L by b becomes 0 and a partial differential of L by a becomes 0, following relations can be obtained.

\begin{align} &\cfrac{\partial L}{\partial b}=0 & \rightarrow & a=\cfrac{\sum_{j=1}^N \sqrt{b x_j} - 2 N}{\sum_{j=1}^N \left[(b x_j)\cdot \exp\left(-\sqrt{b x_j}\right)\right]} &= a_1 \\ &\cfrac{\partial L}{\partial a}=0 & \rightarrow & a=\cfrac{N}{\sum_{j=1}^N \exp\left(-\sqrt{b x_j}\right) + \sum_{j=1}^N \left[\sqrt{b x_j}\cdot \exp\left(-\sqrt{b x_j}\right)\right]} &= a_2 \end{align}

Since a function L has the maximum value at the condition of a_1=a_2, the parameter b can be obtained using Bisection Method considering following equation.

\begin{equation} h(b)=a_1(b)-a_2(b) \end{equation}

It shall be noted that a relation of a_2 \gt 0 is everytime true, but b has to satisfy the following relation to ensure the condition of a_1 \gt 0.

\begin{equation} a_1 \gt 0 \quad \rightarrow \quad b \gt \left(\cfrac{2 N}{\sum_{j=1}^N \sqrt{x_j}}\right)^2 \end{equation}

Above relation can be used as the lowest limit of b in Bisection Method. Regarding the highest limit of b for Bisection Method, it is no information. Therefore, following method in which searching range of solution is moved with an interval of 0.5 is applied for the calculation.

# Sample Code for Bisection method by Python
# (xx : vector of observed values)
# (bb : unknown scalar parameter to be solved)
        def func(bb,xx):
            # function for Bisection Method
            nn=len(xx)
            a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
            a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx))))
            return a1-a2
        imax=100    # maximum number of iterations
        err=1.0e-10 # allowable error of the zero value
        x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2
        x2=x1+0.5
        for i in range(0,imax):
            xi=0.5*(x1+x2)
            f1=func(x1,xx)
            f2=func(x2,xx)
            fi=func(xi,xx)
            if 0<f1*f2:
                x1=x2
                x2=x1+0.5
            elif abs(x2-x1) < err:
                break # converged
            else:
                if f1*fi < 0.0: x2=xi
                if fi*f2 < 0.0: x1=xi
        return xi

結果出力

用いたデータは、「設計 確率雨量の推定（３）」で用いた１つである、宮崎の日最大雨量データである。

従来法Scipy 非線形回帰
f:id:damyarou:20190601110417j:plain f:id:damyarou:20190601110507j:plain


Return Period and Maximum Daily Raimfall (unit: mm/day)
Year 2 5 10 20 50 100 200 500 1000 2000 5000 10000
Conventional 166.7 230.3 277.4 326.2 394.8 450.0 508.4 590.4 656.1 725.1 821.2 897.6
Scipy leasrsq 166.7 237.5 290.4 345.6 423.5 486.5 553.3 647.4 723.1 802.6 913.71002.2

従来法によるプログラム

上述した、二分法により未知パラメータ a および b を計算する方法によるプログラムを以下に示す。 プロッティング・ポジション公式は、Hazen 公式としている。

import numpy as np
import matplotlib.pyplot as plt


def qqplot(xx,ye,ss):
    fsz=12
    xmin=0; xmax=700; dx=100
    ymin=0; ymax=700; dy=100
    plt.figure(figsize=(4,4),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Observed (mm/day)')
    plt.ylabel('Predicted (mm/day)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.plot(xx,ye,'o',clip_on=False)
    plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1)
    xs=xmin+0.05*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished')
    plt.tight_layout()
    fnameF='fig_qq_M_org.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()    

    
def sqrtet(xx,pp,pyear):
    def est_sqrtet(aa,bb,pp,x0):
        # estimation of quantile by Newton-Raphson Method
        def gx(x,aa,p):
            ggx=np.log(1+x)-x-np.log(-1/aa*np.log(p))
            gdx=1/(1+x)-1
            return ggx,gdx
        imax=100
        eps=1.0e-10
        tp=np.zeros(len(pp),dtype=np.float64)
        x1=x0
        for i,p in enumerate(pp):
            for j in range(imax):
                ggx,gdx=gx(x1,aa,p)
                x2=x1-ggx/gdx
                if np.abs(x2-x1)<eps: break
                x1=x2
            tp[i]=x2
        return tp**2/bb
    
    def bisection(xx):
        # estimation of coefficient bb by Bisection Method
        def func(bb,xx):
            # function for Bisection Method
            nn=len(xx)
            a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
            a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx))))
            return a1-a2
        imax=100    # maximum number of iterations
        err=1.0e-10 # allowable error of the zero value
        x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2
        x2=x1+0.5
        for i in range(0,imax):
            xi=0.5*(x1+x2)
            f1=func(x1,xx)
            f2=func(x2,xx)
            fi=func(xi,xx)
            if 0<f1*f2:
                x1=x2
                x2=x1+0.5
            elif abs(x2-x1) < err:
                break # converged
            else:
                if f1*fi < 0.0: x2=xi
                if fi*f2 < 0.0: x1=xi
        return xi    
    
    bb=bisection(xx)
    aa=(np.sum(np.sqrt(bb*xx))-2*len(xx))/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
    # estimation
    x0=np.sqrt(bb*xx[-1])
    ye=est_sqrtet(aa,bb,pp,x0)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_sqrtet(aa,bb,prob,x0)
    res=np.array([aa,bb,rr])
    ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss


def main():
    pyear=np.array([2,5,10,20,50,100,200,500,1000,2000,5000,10000])
    xx=np.array([490.2, 175.7, 285.6, 152.9, 199.1,  98.4, 135.3, 188.3, 143.1, 157.2,
                 165.1, 119.5, 175.9, 198.7, 293.0, 122.4, 104.3, 154.8, 112.0, 305.8,
                 108.6, 178.8, 135.3, 136.5, 142.9, 175.4, 197.2, 165.0, 232.5,  97.4,
                 328.9, 316.7, 125.1, 194.9, 232.2, 182.8, 298.1, 115.9, 175.4, 180.1,
                  99.8, 321.7, 146.7, 156.7, 191.5, 163.2,  83.6, 217.9, 430.4, 143.7,
                 137.6, 132.8, 202.9, 587.2, 157.5, 170.1, 213.4, 301.9, 112.9, 205.2,
                 168.0, 129.5, 175.8, 155.6, 221.4, 341.9, 171.5, 181.1, 173.0, 151.5,
                 116.3, 145.4, 124.6, 172.7, 135.0, 197.0, 124.8, 172.6, 207.7, 131.1,
                 387.3,  98.9, 242.0, 112.5, 100.0, 200.0, 186.5, 108.0, 132.5, 162.5,
                 121.0, 155.5,  93.5, 286.5, 167.0, 103.5, 141.5, 315.5, 102.5, 207.5,
                  86.0, 233.5, 146.5, 228.5, 437.5, 108.5, 120.0, 317.5, 139.5, 203.0,
                 132.5, 234.0, 144.5, 217.5, 133.0, 217.0,  91.0, 135.5, 188.0, 310.0,
                 155.0, 199.5, 280.5, 137.5, 125.5, 150.0, 247.5, 119.0, 200.0, 124.5,
                 222.0, 268.0]) # Miyazaki

    # coefficient for plotting position formula
    #alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    alpha=0.5 # Hazen
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    res,ye,yp,ss=sqrtet(xx,pp,pyear)    
    qqplot(xx,ye,ss)

    s0='{0:4s}'.format('Year')
    s1='{0:4s}'.format('SQRT')
    for i in range(len(pyear)):
        s0=s0+',{0:6.0f}'.format(pyear[i])
        s1=s1+',{0:6.1f}'.format(yp[i])
    print(s0)
    print(s1)
        
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Scipyの非線形最小二乗法による回帰プログラム

未知パラメータ a および b を、Scipy の非線形最小二乗法による回帰計算で求めるプログラムを示す。 ここで、非線形回帰を行うために a および b の初期値を入力する必要があるが、ここでは以下のように初期値を定めた。

b の初期値

b の初期値としては小さめの値を与えたいため、二分法の解の検索区間の最小値である以下の値を b の初期値とした。

\begin{equation} b_0=\left(\cfrac{2 N}{\sum_{j=1}^N \sqrt{x_j}}\right)^2 \end{equation}

a の初期値

クオンタイル x_p を求めるための t と p および a の関係式より、以下の関係を導ける。

\begin{equation} \ln(1+t_p) - t_p = \ln\left[ -\cfrac{1}{a}\ln(p) \right] \quad \rightarrow \quad a = -\cfrac{\ln(p)}{\exp[\ln(1+t)-t]} \end{equation}

上式において、大き目の a を初期値として与えたいため、 p=1\times 10^{-6}、 t=\sqrt{b_0 \cdot x_{max}} としたときの a を初期値とした。すなわち、

\begin{equation} a_0 = -\cfrac{\ln(1\times 10^{-6})}{\exp[\ln(1+t_0)-t_0]} \qquad t_0=\sqrt{b_0 \cdot x_{max}} \end{equation}

プログラム全文

プログラム全文を以下に示す。プロッティング・ポジション公式は、Hazen 公式としている。

import numpy as np
from scipy import optimize # for use of leastsq
import matplotlib.pyplot as plt


def qqplot(xx,ye,ss):
    fsz=12
    xmin=0; xmax=700; dx=100
    ymin=0; ymax=700; dy=100
    plt.figure(figsize=(4,4),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Observed (mm/day)')
    plt.ylabel('Predicted (mm/day)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.plot(xx,ye,'o',clip_on=False)
    plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1)
    xs=xmin+0.05*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished')
    plt.tight_layout()
    fnameF='fig_qq_M_reg.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()    

    
def sqrtet(xx,pp,pyear):
    def est_sqrtet(aa,bb,pp,x0):
        # estimation of quantile by Newton-Raphson Method
        def fx(tp,aa,p):
            return np.log(1+tp)-tp-np.log(-1/aa*np.log(p))
        tp=np.zeros(len(pp),dtype=np.float64)
        for i,p in enumerate(pp):
            root = optimize.newton(fx,x0,args=(aa,p))
            tp[i]=root
        return tp**2/bb

    def func_sqrtet(parameter,xx,pp):
        # function for leastsq
        aa=parameter[0]
        bb=parameter[1]
        x0=np.sqrt(bb*xx[-1])
        est=est_sqrtet(aa,bb,pp,x0)
        residual=xx-est
        return residual
    
    # parameter calculation by scipy.optimize.leastsq
    b0=(2*len(xx)/np.sum(np.sqrt(xx)))**2
    t0=np.sqrt(b0*xx[-1])
    a0=-np.log(1e-6)/np.exp(np.log(1+t0)-t0)
    npara=2; parameter0 = np.array([a0,b0])
    result = optimize.leastsq(func_sqrtet,parameter0,args=(xx,pp))
    aa=result[0][0]
    bb=result[0][1]
    # estimation
    x0=np.sqrt(bb*xx[-1])
    ye=est_sqrtet(aa,bb,pp,x0)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_sqrtet(aa,bb,prob,x0)
    res=np.array([aa,bb,rr])
    ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss


def main():
    pyear=np.array([2,5,10,20,50,100,200,500,1000,2000,5000,10000])
    xx=np.array([490.2, 175.7, 285.6, 152.9, 199.1,  98.4, 135.3, 188.3, 143.1, 157.2,
                 165.1, 119.5, 175.9, 198.7, 293.0, 122.4, 104.3, 154.8, 112.0, 305.8,
                 108.6, 178.8, 135.3, 136.5, 142.9, 175.4, 197.2, 165.0, 232.5,  97.4,
                 328.9, 316.7, 125.1, 194.9, 232.2, 182.8, 298.1, 115.9, 175.4, 180.1,
                  99.8, 321.7, 146.7, 156.7, 191.5, 163.2,  83.6, 217.9, 430.4, 143.7,
                 137.6, 132.8, 202.9, 587.2, 157.5, 170.1, 213.4, 301.9, 112.9, 205.2,
                 168.0, 129.5, 175.8, 155.6, 221.4, 341.9, 171.5, 181.1, 173.0, 151.5,
                 116.3, 145.4, 124.6, 172.7, 135.0, 197.0, 124.8, 172.6, 207.7, 131.1,
                 387.3,  98.9, 242.0, 112.5, 100.0, 200.0, 186.5, 108.0, 132.5, 162.5,
                 121.0, 155.5,  93.5, 286.5, 167.0, 103.5, 141.5, 315.5, 102.5, 207.5,
                  86.0, 233.5, 146.5, 228.5, 437.5, 108.5, 120.0, 317.5, 139.5, 203.0,
                 132.5, 234.0, 144.5, 217.5, 133.0, 217.0,  91.0, 135.5, 188.0, 310.0,
                 155.0, 199.5, 280.5, 137.5, 125.5, 150.0, 247.5, 119.0, 200.0, 124.5,
                 222.0, 268.0]) # Miyazaki

    # coefficient for plotting position formula
    #alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    alpha=0.5 # Hazen
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    res,ye,yp,ss=sqrtet(xx,pp,pyear)    
    qqplot(xx,ye,ss)

    s0='{0:4s}'.format('Year')
    s1='{0:4s}'.format('SQRT')
    for i in range(len(pyear)):
        s0=s0+',{0:6.0f}'.format(pyear[i])
        s1=s1+',{0:6.1f}'.format(yp[i])
    print(s0)
    print(s1)
        
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く