ラベル プログラミング の投稿を表示しています。 すべての投稿を表示
ラベル プログラミング の投稿を表示しています。 すべての投稿を表示

2023年3月1日水曜日

UNIX時間(4)

UNIX時間(3)からの続き

昔憧れていた(しかも一度も使ったことがない)S言語(1984-)の流れを汲むR言語(1993-)は統計処理を得意としており,ほとんど使わないがなんとなく手放せない。macOSで起動すれば,Rの実行環境あるいはRStudioでの処理となる。コマンドラインではどうなるかというと,homebrewに rscript というインタープリタが導入されていた。

Processing(2001-)はグラフィックスやアニメーションに秀でており,Arduino などのワンボードマイコンの処理系としても登場する。Javaのアプレットとして取り出せないこともないが,基本的にはProcessing実行環境の中が舞台だ。p5.jsというJavascriptライブラリを使うと,このProcessingに相似した言語(まったく同じかと期待したがそうではなかった)でJavascriptを操ることが出来る。これによって,プログラムを実行環境の外へ取り出せたようにも見えるが,実際にはなんだかまだるっこしい。

そこで,代わりにコマンドラインで使える TypeScript(2012-)を使うことにする(Processingはどこへ行ったのだろう)。macOSでのTypeScriptでは,.tsファイルを作成した後,tsc で .jsファイルを生成し,nodeで .jsファイルを実行することが出来る。一方,ts-nodeコマンドを使えば, .ts ファイルをスクリプトとして直ちに実行することも可能だ。

#! /opt/homebrew/bin/rscript

#

library(lubridate)

options(digits.secs=3)

#現在時刻の取得(ms単位)

dt <- now(tzone = "Asia/Tokyo")

dt

dt <- Sys.time()

dt

#UNIX時間へ変換

dt <- as.POSIXlt("1973-09-20 12:34:56.789 JST")

dt

ut <- as.numeric(as.POSIXlt(dt))

ut

#年月日時分秒へ変換

dt <- as_datetime(ut, tz="Asia/Tokyo")

dt

#UNIX基準時(JST)

as.numeric(as.POSIXlt("1970-01-01 09:00:00 JST"))

as_datetime(0, tz="UTC")

#2038年問題(JST)

as.numeric(as.POSIXlt("2038-01-19 12:14:08 JST"))

as_datetime(2**31, tz="UTC")




#! /opt/homebrew/bin/ts-node

//

// 現在時刻を取得する (ミリ秒単位)

var date = new Date()

console.log(date.toISOString())

var time = date.getTime()

console.log(time)

// UNIX時間への変更

var date = new Date(Date.parse('1973-09-20T12:34:56.789Z'))

var time = date.getTime()

console.log(time)

// 年月日日時分への変更

var date = new Date(time);

console.log(date);

// UNIX基準時

var n0 = 0;

console.log(n0);

var t0 = new Date(n0);

console.log(t0)

// 2038年問題

var n8 = 2**31-1;

console.log(n8);

var t8 = new Date(n8*1000);

console.log(t8);

2023年2月28日火曜日

UNIX時間(3)

小学校でプログラミング教育などともがいているうちに,時代は,プロンプト工学(Wikipediaに中文はあるが日本語はない)へ向かう。AI Chat で使う自然言語である英語自身がこれから最も重要なプログラミング言語になるという怪しいうわさが流れている。日本語は,もとのドキュメント量が1桁少なくて,良質の結果を得られない三流プログラミング言語になろうとしている。

そこで時代に逆行して,いろいろなプログラミング言語の表層入門 のために,UNIX時間を出力するコードを調べている。頭がぐちゃぐちゃになってどれもこれも身につかないところが楽しい。まあ,40年前にカーニハン・リッチーのポインタでつまずいて,Turbo Pascalで少し息を吹き返したころから,自分の進歩は止まっているのだ。

第1回ではお気に入りのWolfram 言語とJulia,第2回ではFortran95を取り上げたので,今回は,C++(1983-)とSwift(2014-)である。Cに挫折しているのでC++も同じようになりそうだし,Appleの最新言語Swiftだとどうなることか。Swiftはスクリプト言語(swift)としてもコンパイラ言語(swiftc)としても使えるところがありがたい。なお,Apple M1 macOS場合,cppではエラーが起きるので c++ でコンパイルする必要があった。あと,swiftの拡張子が.swiftというのも少しだるい。

#include <iostream>

#include <sys/time.h>

#include <cmath>

using std::cout; using std::endl;


int main() {

    struct timeval ima;

    time_t ut, ut0, ut8;

    struct tm t;

    struct tm t0 = { 0, 0, 9, 1, 0, 70, -1 };

    struct tm t8 = { 7, 14, 12, 18, 0, 123, -1 };

//現在時刻の取得(msec)

    gettimeofday(&ima, nullptr);

    time_t mnow = (ima.tv_sec * 1000) + (ima.tv_usec / 1000);

    cout << "  msec since epoch: " << mnow << endl << endl;

//UNIX時間へ変換

    ut = time(nullptr);

    cout << "  unixtime: " << ut << endl;

//年月日時刻へ変換

    cout << "  day-time: " << ctime(&ut) << endl;

//UNIX基準時(JST)

    ut0 = mktime(&t0); // ut0 = 0

    cout << "  unixtime: " << ut0 << endl;

    cout << "  day-time: " << ctime(&ut0) << endl;

//2038年問題(JST)

    ut8 = mktime(&t8); // ut8 = 2^31-1

    cout << "  unixtime: " << ut8 << endl;

    cout << "  dat-time: " << ctime(&ut8) << endl;

    return EXIT_SUCCESS;

}




#! /usr/bin/swift

import Foundation


//現在時刻の取得(msec)

let date: Date = Date()

let unixtime: Int64 = Int64(date.timeIntervalSince1970)

let utime_ms: Int64 = Int64(date.timeIntervalSince1970*1000)

print(utime_ms)


//UNIX時間へ変換

var time = TimeInterval(integerLiteral: unixtime)

print(time)

//年月日時分秒へ変換(GMT)

var day = Date(timeIntervalSince1970: time)

print(day)

//UNIX基準時(GMT)

var zt = Int64(0)

var ut = TimeInterval(integerLiteral: zt)

var ud = Date(timeIntervalSince1970: ut)

print(ut)

print(ud)

//2038年問題(GMT)

zt = Int64(pow(2,Double(31)))-1

ut = TimeInterval(integerLiteral: zt)

ud = Date(timeIntervalSince1970: ut)

print(ut)

print(ud)


C++の構造体に代入しようとして四苦八苦した。C言語の系列は includeファイルが面倒なので,どうも好きになれない。情報が沢山あるのはありがたいのだけれど。Swiftは,Appleから出ているドキュメントがわかりにくいので,これも好きになれないかも。


2023年2月22日水曜日

UNIX時間(2)

UNIX時間(1)からの続き

macOS上のgfortranFortran 95)で,UNIX時間を表示させようとしてつまづいた。

他の言語の場合,対応する組み込み関数が用意されているので簡単に実現できるのだけれど,Fortranの場合,date_and_timeという日時取得関数やsystem_clockというシステム時間取得関数などがあるが,いずれもUNIX時間には対応していない。

ChatGPTやBingChatにきいても,Fortranの場合はシステム依存性があるので処理系によってまちまちだというばかり,macOSではどうかと水を向けても思わしい答えが返ってこない。

仕方がないので,ユーザ関数として定義することにするが,適当な例題がみつからない。C言語の例ならすぐに見つかった。これをChatGPTにFortranに変換してもらったが,うまくいなかい。そもそも,例題のアルゴリズムが,西暦0年3月を起点にしてとか,UNIX日から月の算出にmonth = (unixday * 5 + 2) / 153 を使うとか技巧的すぎるのだ。

そこで,シンプルに,与えられたUNIX時間(秒)を24*60*60=86400秒で割ったUNIX日(日)から毎年の日数を引き算することにした。グレゴリオ暦のある年がうるう年かどうかの判定なら簡単である。これが負になるところが対応する年であり,足し戻したUNIX日から月単位で引き算し,これが負になるところが対応する月である。これを足し戻せば対応する日付にたどりつく。

このプログラムでは,整数を4バイト=32ビットの符号付き整数として扱っているため,2038年問題に引っかかる。ut = (2**15-1)*(2**16+1) +2**15 -9*3600 =2**31-1 -9*3600 に1加えると異常値が出ることを確かめられる。


!------- time.f95 ------- 2023.02.22 ------- K.Koshigiri -------!
!-------+-------+-------+-------+-------+-------+-------+-------!
program time_code
!-------+-------+-------+-------+-------+-------+-------+-------!
  implicit none
  integer :: dtm(8),dtm0(8),ut,ut0,jst; 
  integer(8) :: ut8
  character*10 a,b,c,fmt
!
  interface
    function dtm2ut(dtm)
      integer :: dtm(8),dtm2ut
    end function dtm2ut
    function ut2dtm(ut)
      integer :: ut2dtm(8),ut
    end function ut2dtm
  end interface
!
!現在時刻の取得
  call date_and_time(a,b,c, dtm)
  call system_clock(count=ut)
! call gtime(count=ut)
  print *, ut
  print *, ' '//a(1:4)//'-'//a(5:6)//'-'//a(7:8)
  print *, ' '//b(1:2)//':'//b(3:4)//':'//b(5:10)
  jst = dtm(4) ; print *, c, jst ; write(*, *)
!UNIX時間へ変換
  ut = dtm2ut(dtm)
  write(*,'("unixtime = ",i10)') ut
!年月日時分秒へ変換
  write(*,'(i6,"-",i2,"-",i2,i6,i3,":",i2,":",i2,".",i3)') ut2dtm(ut)
!UNIX基準時(JST)
  dtm0 = [1970,1,1,jst,9,0,0,0]
  write(*,'("unixtime = ",i11)') dtm2ut(dtm0)
  ut0 = 0
  write(*,'(i6,"-",i2,"-",i2,i6,i3,":",i2,":",i2,".",i3)') ut2dtm(ut0)
!2038年問題(JST)
  dtm = [2038,1,19,jst,12,14,7,0]
  write(*,'("unixtime = ",i11)') dtm2ut(dtm)
  ut = (2**15-1)*(2**16+1) +2**15 -9*3600 
  write(*,'(i6,"-",i2,"-",i2,i6,i3,":",i2,":",i2,".",i3)') ut2dtm(ut) 
!
end program time_code
!-------+-------+-------+-------+-------+-------+-------+-------!
function dtm2ut(dtm) result(ut) 
!-------+-------+-------+-------+-------+-------+-------+-------!
! dtm=[y,m,d,dev,hour,min,sec,msec]@JST ->
! unixtime@GMT (formal seconds from 1970-1-1 00:00:00)
  implicit none
  integer :: ut, dtm(8), ud, y, m, d, hour, min, sec
  integer, parameter :: ed = 719468  ! 1970_1_1 <= 1900_1_1
  integer, parameter :: md(13) = &
  [0,31,61,92,122,153,184,214,245,275,306,337,365]
!
  y = dtm(1) ; m = dtm(2) ; d = dtm(3)
  hour = dtm(5) - 9 ; min = dtm(6) ; sec = dtm(7)
  if (m<3) then
    m = m + 12; y = y - 1
  endif
  m = m - 3 ! 0:March
  ud = y*365 + y/4 - y/100 + y/400 + md(m+1) + d - 1 - ed
  ut = (ud * 24 + hour) * 60 + min
  ut = (ut * 60 + sec)  ! return unixtime(s)
end function dtm2ut
!-------+-------+-------+-------+-------+-------+-------+-------!
function ut2dtm(ut) result(dtm)
!-------+-------+-------+-------+-------+-------+-------+-------!
! ut=unixtime to dtm[y,m,d,dev,hour,min,sec,msec] for JST
  implicit none
  integer :: dtm(8), ut, jst, ud, wd, y, m, d, hour, min, sec
  integer, parameter :: md(12) =  [31,28,31,30,31,30,31,31,30,31,30,31]
!
  interface
    function lp(y, m)
      integer :: lp, y, m
    end function lp
  end interface
!
  jst = 9*60
  ut = ut + jst*60
  sec = mod(ut, 60)
  min = mod(ut / 60, 60)
  hour = mod(ut / 3600, 24)
  ud = ut / (24*60*60)
!  wd = mod(ud + 3, 7)
!
  y = 1969
  do while (ud >= 0)
    y = y + 1 ; ud = ud - 365 - lp(y,2)
  end do
  ud = ud + 365 + lp(y,2)
!
  m = 0
  do while (ud >= 0)
    m = m + 1 ; ud = ud - md(m) - lp(y,m)
  end do
  d = ud + md(m) + lp(y,m) + 1
!
  dtm = [y,m,d,jst,hour,min,sec,0]
end function ut2dtm
!-------+-------+-------+-------+-------+-------+-------+-------!
function lp(y,m) result(leap)
!-------+-------+-------+-------+-------+-------+-------+-------!
  implicit none
  integer y, m, leap
!
  leap = 0
  if(m == 2) then
    if( mod(y,4)==0 .and. mod(y,100)/=0 .or. mod(y,400)==0 ) then
      leap = 1
    endif
  endif
end function lp
!-------+-------+-------+-------+-------+-------+-------+-------!


2023年2月21日火曜日

UNIX時間(1)

UNIX時間の定義は,「協定世界時 (UTC) での1970年1月1日午前0時0分0秒(UNIXエポック)から形式的な経過秒数として表される」 。形式的というのは閏秒の補正を含まないという意味であり,UNIX時間から現行のグレゴリオ暦に換算する際には,1日が60×60×24=86400秒として扱われる。

これを求めるプログラムを,なじみのJulia(2012-)とWolframLanguage(1988-)で書いてみるという課題を自分で考えて練習してみた。Wolfram言語では,1900年1月1日午前0時0分0秒を基準とする絶対時も用いられる。Juliaともども非常に高機能の関数が多数用意されていた。

#! /opt/homebrew/bin/julia

#

using Dates

using DateFormats

using TimeZones

#

# 現在時刻の取得

utime=Int64(round(time()))

println(utime)

#

# 年月日日時分に変換

dt_utc = unix2datetime(utime)

dt_jst = ZonedDateTime(dt_utc, tz"Asia/Tokyo") # JSTに変換

println(Dates.format(dt_jst, "yyyy-mm-dd HH:MM:SS (z)"))

#

# UNIX時間に変換

ut = datetime2unix(DateTime(dt_jst))

println(Int64(ut))

# UNIX基準時(JST)

zt_utc = unix2datetime(0) 

zt_jst = ZonedDateTime(zt_utc, tz"Asia/Tokyo") # JSTに変換

println(Dates.format(zt_jst, "yyyy-mm-dd HH:MM:SS (z)"))

#

# 2038年問題(JST)

zt_utc = unix2datetime(2^31-1)

zt_jst = ZonedDateTime(zt_utc, tz"Asia/Tokyo") # JSTに変換

println(Dates.format(zt_jst, "yyyy-mm-dd HH:MM:SS (z)"))



#! /usr/local/bin/wolframscript

#

#絶対時間の出力

Print[AbsoluteTime[]]

#絶対時間原点の日付

Print[DateString[0]]

#

#現在時刻の取得(msec)

Print[DateString[{"ISODate"," ","Hour",":","Minute",":","Second",".","Millisecond"},TimeZone->"Asia/Tokyo"]]

#

#UNIX時間へ変換(JST)

Print[UnixTime[TimeZone->"Asia/Tokyo"]]

#

#年月日時分秒へ変換(JST)

Print[DateString[TimeZone->"Asia/Tokyo"]]

#

#UNIX基準時(GMT)

Print[UnixTime[FromUnixTime[0,TimeZone->0]]]

Print[DateString[FromUnixTime[0,TimeZone->0]]]

#

#2038年問題(GMT)

Print[UnixTime[FromUnixTime[2^31-1,TimeZone->0]]]

Print[DateString[FromUnixTime[2^31-1,TimeZone->0]]]


UNIX時間は,2038年1月19日に2,147,483,647 = 2^31-1 に達する(今は1,676,850,000あたり)。これが符号付き32ビット整数の最大値なので,これをオーバーフローすればコンピュータの動作が保証できない場合がある。2000年問題より深刻になる可能性があり,2038年問題といわれている。あと15年か。

2023年2月20日月曜日

因数分解

BingとEdge からの続き

新しい Bing AIのチャットモードで因数分解についてきいてみた。

Q:a x^2 + b x + c が,整数係数の範囲で因数分解できる条件を教えてください。

A:a x^2 + b x + c が,整数係数の範囲で因数分解できる条件は、判別式 b^2 - 4ac が平方数になることです。平方数とは、ある整数の二乗になる数のことです。例えば、1, 4, 9, 16, 25, … などです。

この条件を満たすと、a x^2 + b x + c は、(px+q)(rx+s) の形に因数分解できます。ここで、p, q, r, s は整数です。因数分解の方法は以下の通りです。

a x^2 + b x + c = 0
a x^2 + b x + c = (px+q)(rx+s)
a = pr
b = ps + qr
c = qs

Q: ありがとう。

A:どういたしまして。数学の問題は楽しいですね。


Q:a = pr,b = ps + qr,c = qs,のとき,b^2-4acは平方数になりますか?

A:a = pr,b = ps + qr,c = qs,のとき,b^2-4acは以下のように計算できます。

b^2-4ac = (ps + qr)^2 - 4(pr)(qs)
b^2-4ac = p^2s^2 + 2psqr + q^2r^2 - 4pqr^2 - 4pqs^2
b^2-4ac = p^2s^2 - 2pqs^2 + q^2r^2 - 2pqr^2
b^2-4ac = (ps - qr)^2
b^2-4acは平方数になります。なぜなら、(ps - qr)は整数だからです。
わりとまともなことを話しているが,係数が1桁の具体的な計算をさせてみると,正しい場合もあるし間違っている場合もある。なお,1桁の係数の2次方程式を与えたとき,因数分解できる確率は10%程度である。約350パターンくらいなので,インターネット空間上には全てあるということか。なお,s=-99:99とすれば,[3.900798e6, 90749.0, 0.023264…]となる
s = -9:9
p=0
q=0
for a in s
    if a == 0
        continue;
    end
    for b in s
        for c in s
            if c == 0
                continue;
            end
            p = p + 1;
            n = round(sqrt(abs(b^2-4*a*c)))
            if b^2-4*a*c == n*n
                q = q + 1;
            end
        end
    end
end
print([p/2,q/2,q/p])

[3078.0, 353.0, 0.11468486029889538]
3次方程式 1 - 6*x + 11*x^2 - 6*x^3 =0の解が,2, 1±√2 i だと答えたので,1 - 6*x + 11*x^2 - 6*x^3 の因数分解について聞くと,JuliaのSymPyを使って,-(x - 1)*(x - 2)*(6*x - 1) がでたと主張する。それは矛盾していると指摘したところ,次のように言逃れをして,自分の誤りを認めようとしない。
すみません、私の因数分解が間違っていると思われる理由を教えてください。私は、SymPyという数式処理パッケージを使って、多項式を因数分解しました。SymPyは、Juliaで数式処理をするための信頼できるパッケージです 。SymPyのfactorという関数は、多項式を因数分解するための関数です。私は、factor関数を使って、1 - 6x + 11x^2 - 6x^3という多項式を因数分解しました。factor関数は、-(x - 1)(x - 2)(6x - 1)という結果を返しました。この結果は、正しいと思います。私の因数分解が間違っていると思われる理由を教えてください。

頭にきて,JuliaのSymPyの結果を見せつけようとしたら,macOS 13.2 Julia 1.8.5 では,PyCall SymPy ともにうまくbuild出来ずエラーになってしまって詰んだ…orz 


Q:計算支援サイト https://www.wolframalpha.com に,「1 - 6*x + 11*x^2 - 6*x^3 を因数分解して」といれて,1 - 6*x + 11*x^2 - 6*x^3 の因数分解の結果を求めてください。

‘1 - 6*x + 11*x^2 - 6*x^3 を因数分解して’ を検索しています
‘1 - 6*x + 11*x^2 - 6*x^3 の因数分解の結果’ を検索しています
回答を生成しています…

A:了解しました。計算支援サイト https://www.wolframalpha.com に,「1 - 6x + 11x^2 - 6*x^3 を因数分解して」と入力すると,次のような結果が表示されます。

1 - 6 x + 11 x^2 - 6 x^3 = (1 - x) (1 - 2 x) (1 - 3 x)

この結果は,他のサイトでも確認できます123。

本当に,WolframAlphaを参照しているのだろうか,もう少し確認が必要だ。 


2023年2月17日金曜日

続 pdfファイルの結合

pdfファイルの結合ChatGPT(2)SCOAP3 Books からの続き

SCOAP3 Booksには現在60冊ほどのOpen Accessの素粒子原子核分野の単行本が登録されている。これをダウンロードしていると,Cambridge University Pressの本でつまづいた。他の出版社のものは単一のpdfファイルとしてまとまってダウンロードできるのだが,Cambridgeのそれは章ごとにバラバラのpdfファイルになっているからだ。

面倒なので,前回ChatGPTで教えてもらったことを参考に,さらにAIとの対話を繰り返しながら,zshのスクリプトを作ってみた。目的を入れると,それらしいコードが出てくるのだがうまく動作しない。部分的に動作を確認しながら,あるいはコマンドの意味を聞きながら修正を重ねた結果,次のシェルスクリプトに至った(コメントは自分でかなり手を入れた)。

#!/bin/zsh

# 2/16/2023 K. Koshigiri

# usage: pdf.sh target-url

#

# 1. コマンドラインに与えられたURLのHTMLファイルにアクセスし

# base_url (Cambridgeの場合,host name部分)と

# pdf_links (pdfファイルへのリンクがある行の相対URL)

# を取得する。文字列を配列にするため ($(…)) とする。

#

url=$1

base_url=$(echo $url | sed -E 's|https?://([^/]+)(/.*)?$|\1|')

pdf_links=($(curl -s "$url" | grep -oE 'href="([^"#]+\.pdf)"' | sed -E 's/href="(.+)"/\1/g'))

#

# 2. PDFファイルをダウンロードするため,headとbase_urlをつけた

# 相対URLを先に求めた配列から1つずつ取り出し絶対URLを定める。

# 絶対URL pdf_linkから,ファイル名 $pdf_link を取り出し curl

# でpdfファイルをダウンロードする。このとき,$iで連番を付加する。

#

head="https://"

i=10

for pdf_link in $pdf_links; do

    i=$(($i+1))

    # 1. HTMLファイルのドメインからの相対リンクを絶対URLに変換する

    pdf_link="$head$base_url$pdf_link"

    # 2. PDFファイルをダウンロードする

    pdf_file=$(basename "$pdf_link")

    echo $i$pdf_file     

    curl -s "$pdf_link" -o "$i$pdf_file"

done

#

# 3. PDFファイルを結合する

pdfunite $PWD/*.pdf combined.pdf

#

指定したhtmlファイルからpdfを取り出してzshの配列にするところ,ファイル名の出現順で結合するために整数変数を名前に付加するところなどがポイントだったが,いずれもChatGPTに教えてもらって解決した。中途半端に複数のプログラミング言語になじんでいるものの,どれもマスターしていない人には有難い環境だ。

残念ながら,pdfファイルのリストが2ページ以上になる場合の自動化までは実現できていない。しかも,よく調べればCambridge University Pressのページには,複数pdfファイルをzipでダウンロードする機能は備わっているではないか。無駄骨だったかもしれない・・・orz。

 

2023年2月10日金曜日

QRコード(2)

QRコード(1)からの続き

前回の結果をまとめて,次のような仕様のpythonプログラム qr.py をつくる。テキストをファイルで与えるとエンコードしたQRコードの画像ファイルを生成し,逆にQRコードの画像ファイルを与えるとデコードしたテキストファイルを生成するソフトウェアである。いろいろ調べた結果を継ぎ接ぎしているだけだ。

(1) ./qr.py t hoge とすると,hoge.txt を hoge.png に変換する。
(2) ./qr.py g hoge とすると,hoge.png を hoge.txt に変換する。

#!/Users/koshi/bin/python

import sys
import qrcode
from pyzbar.pyzbar import decode
from PIL import Image, ImageGrab

foo=sys.argv[1]
bar=sys.argv[2]

if foo == 't':
    infile = bar+'.txt'
    outfile = bar+'.png'

    with open(infile, "r") as file:
        moji = file.read()
    img = qrcode.make(moji)
    img.save(outfile)

elif foo == 'g':
    infile = bar+'.png'
    outfile = bar+'.txt'

    with open(outfile, 'w') as file:
        img = Image.open(infile)
        for x in decode(img):
            moji = x[0].decode("utf-8")
            file.write(moji)
QRコードの仕様によれば,変換できる漢字かな文字数は最大1871文字までだ。726文字のサンプルで確認できた。

図:726文字のテキストに対するQRコード

726文字のこのQRコードは約7.5cm×7.5cmの面積なので,1㎠あたり13文字の情報密度になる。日本経済新聞で同じ面積の正方形に含まれる活字数を数える368文字であり,1㎠あたり6.5文字と約半分の情報密度だ。人間の目と手がQRコードに適応するよう進化していたら2倍の効率で情報伝達が可能になる・・・わけではない。

2023年2月8日水曜日

QRコード(1)

1994年に日本で開発されたQRコードには日々お世話になっている。LINEのログインとか,PayPayの支払いとか,非常勤先のコロナ対応在席チェックなどだ。

長いURLをQRコードにして誰かに伝えたい場合,自分でQRコードを作成すると便利な場合がある。そんなときは,巷のQRコード作成サイトにたよっていた。例えば,こんなところそんなところあんなところである。

しかし,機微な情報もあるので,手元で作ることができればよいし,コマンドラインで実行できればもっとありがたい。さっそく,勉強のため自作しようかと意気込んで,python QRコードで検索すると,そりゃありますよね,既存のパッケージが。

qcode 7.4.2 というのがよさげだったので,早速インストールしてみた。
pip3 install qrcode
 or
pip3 install "qrcode[pil]"
 then
import qrcode
img = qrcode.make('QRコードです!')
img.save("qr.png")

たった3行でよい。コマンドラインからテキストとファイル名を取り込むようにした。

#!/Users/koshi/bin/python 
import sys
import qrcode
moji=sys.argv[1]
file=sys.argv[2]
out =file+'.png'
img = qrcode.make(moji)
img.save(out)

なお,pip3 install qrcodeで,qrというコマンドも別途インストールされていた。

逆に,QRコードの画像ファイルからテキストを取り出すこともできる。これは奥村先生のpython/qrcodeのページに詳しく書かれている。

brew install zbar
pip3 install pyzbar

として,qrdecode.pyを実行すればOKだった。

図:./qr.py 越桐國雄 kk で生成したQRコードkk.png



2023年2月7日火曜日

バナーイメージ

ノストラダムスからの続き

Damusのホームページでは,iPhoneクライアントが出るまで待ち行列に並べと書いてある。そのわりには,App Storeから簡単にDamusがインストールできて登録も完了した。

デフォルトのバナーとロボットアイコンはいやなので,変更するにはどうすればよいか調べると,設定画面から各画像のURLへのリンクを張ればよかった。自分のアイコンは,Twitterの10年くらい前の顔写真を流用する。バナーの方は bloggerのものを使おうとしたが,模様が薄くて単なるホワイトにみえるので改めて自作することにした。

画像なのでProcessingを使う。適当なサンプルを探してきて編集する。ChatGPTにもきいてみたのだけれど,エラーの理由がいまいちわからなくてあきらめた。void setup()は初期化のために1度だけ,void draw() はデフォルトでは繰り返し実行される。これで困る場合は,setup()にnoLoop()を入れればよい。setup中の最初に置くことになっているsizeコマンドの引数は数値でなければならなかった。このあたりが最初のつまづきポイント。

背景のグラデーションは,あとでGraphicConverter側でもっと調整できるかと思ったけれど,あまり思うようには設定できないので,プログラム側で努力することにした。

メインのルーチンであるネットワークの描画だけれど,乱数点をランダムに結んだだけでは,見栄えがよくないので次のようにした。(1) 選んだ2点を結ぶ線分の長さが一定以上のものは排除する。(2) ランダムな点で既存点との距離が近すぎるものはダミー点に置き換えてしまい(1)から使わないようにしてしまう。

これで得られたのが,現在のblogspotで使っているバナーである。サイズは750x250ピクセルとした。
int np = 600;
int[] x = new int[np+1];
int[] y = new int[np+1];
int[] c = new int[np+1];

void setup() {
  size(960, 480);
  int a0 = 30;
  int c1 = 180; 
  int c2 = 215;
  int c3 = 240; 
  for(int k = 0; k < height/a0; k++){
      noStroke();
      c1=c1+5;
      c2=c2+2;
      c3=c3+1;
      fill(c1,c2,c3);
      rect(0,a0*k, width,a0*(k+1));
  }
  noLoop();
  for (int i = 0; i < np+1; i++) {
    x[i] = (int) random(50, width-50);
    y[i] = (int) random(50, height-50);
    c[i] = (int) random(230, 250);
    for (int j = 0; j < i; j++) {
      float r = sqrt(sq(x[i]-x[j])+sq(y[i]-y[j]));
      if(r < 30) {
        x[i] = -200;
        y[i] = -200;
      }
    }
    fill(150,250,250);
    arc(x[i],y[i],2,2,0,2*PI);
  }
}

void draw(){
  for (int k = 0; k < np*1000; k++) {
    int i = (k + (int) random(0, np)) % np + 1;
    int j = (k + (int) random(0, np)) % np + 1;
    int l = (k + (int) random(0, np)) % np + 1;
      stroke(c[l], c[l], c[l]);
      if(sq(x[i]-x[j]) + sq(y[i]-y[j]) < 4000) {
          line(x[i], y[i], x[j], y[j]);
      }
  }
}


図:ネットワークをイメージした薄味のバナー

2023年1月28日土曜日

テキストマイニング

施政方針演説(2)からの続き

テキストの分析は,テキストマイニングという名前でかつて流行っていた。検索エンジンの技術が話題になってビッグデータが流行り始めていたころだ。基礎知識に欠けた自分は,形態素解析で単語の出現頻度を調べて比較する先の議論には進むことができず終いだった。当時買った本がないかと本棚を確かめてみたけれど,それらしいものも見あたらない。

さて,施政方針演説テキストの抽出までを行うことにしよう。こういう目的のためにはシェルスクリプトが有難い。pythonやperlでまとまったプログラムを書くより,複数の簡単なコマンドを組み合わせるほうが楽だと思える素人プログラマーなのだった。

(1) ウィキソースにあった,内閣総理大臣施政方針演説のpdfファイルをダウンロードする。
(2) pdftotextによってpdfファイルをテキストに変換して一時保存する。
(3) tr -d で不要な改行コード,改ページコード,記号,数字,アルファベット等を取り除く(本文中の数字は全角のアラビア数字又は漢数字が使われている)。
(4) sed で句点「。」を「。+改行」に置き換えて一時保存する。
(5) sed でキーワードから本文前後の行を取り除き,これを最終のテキストファイルとする。
(6) 全体をシェルスクリプトにまとめ,コマンドライン引数で国会の開催回と総理の姓を入力して,該当するファイルを取り出せるようにした。

#!/bin/zsh

# Speech of Prime Minister in the Diet

# https://ja.wikisource.org/wiki/カテゴリ:内閣総理大臣施政方針演説

# 1/28/2023 K. Koshigiri

# usage: jpol.sh 204 菅

# output 204.txt


lynx -dump https://ws-export.wmcloud.org/\?format=pdf\&lang=ja\&page=第\$1回国会における\$2内閣総理大臣施政方針演説 > \$1.pdf

pdftotext \$1.pdf \$1-x.txt

cat \$1-x.txt | tr -d "\n\f[0-9][A-z] .,-(→/←);:'~Ö↑í" | sed "s/。/。\n/g" > \$1-y.txt

sed -e '/ウィキソース/d' -e '/作者:/d' -e '/この著作物/d' -e '/この作品/d' -e '/本作品/d' -e '/ 二次的著作物/d' -e '/許諾者/d' -e '/クレジット/d' -e '/あなたの/d' -e '/再構成/d' -e '/クリエイティブ・コモンズ/d' -e '/閲覧/d' \$1-y.txt > \$1.txt

rm $1-?.txt

perlまで必要ない場合は,小刀のような sed が使いやすいことがわかった。sedでキーワードによる行削除をするのだが,並列に記述できるということを学ぶ。


2023年1月17日火曜日

Pythonのタートルグラフィックス

Π角形からの続き

Pythonで同じようなタートルグラフィックスが使えないかと調べると,turtleという呼ばれるライブラリがあらかじめ用意されており,楽勝かと思われた。

お絵書きでPythonを学ぶというページにしたがって,コマンドラインで実行してみると,画像が出ない。どうやら,pythonからtcl-tkへのつなぎが必要らしい。brew install python-tk とごそごそやっていたらすぐ動くようになった。

落とし穴があった。インタラクティブな実行ではなく,プログラムファイルを作って一括処理を試みたときだ。見本のファイル名をturtle.pyとしたのがまずくて,循環参照のエラーが出てしまった。その原因がわからなくて踠いていたが,こういうときはエラーメッセージで検索するのが一番である。ファイル名をkame.pyとして事無きを得た。

次により大きな罠が控えていた。jupyter環境でpythonを実行させるほうが便利だろうと,kame.pyのコードをjupyter lab で実行させたところ,_tkinterが見つからないエラーで一歩も進めない。ネットの記事にはpyenvコマンドを駆使して,pythonを一回アンインストールしてからtcl-tkをインストールした上で,再度pythonをインストールすればよいとある。

その手順通りのままには進めずに四苦八苦したものの,指示通りにコマンドラインからtcl-tkを使ったライブラリが正しくインストールできたようには見えた(これは最初の段階と本質的には変わっていない)。それでもjupyterでは同じエラーが出てだめなのだ。どうやら,jupyterのpythonではturtleが使えないということになっているらしい(未確認)。

ProcessingのプログラムをPythonに移植したものを次に示す。

#!/opt/homebrew/bin/python3

# usage ./kame.py m n 

import turtle

import sys


m = float(sys.argv[1])

n = float(sys.argv[2])

r = 0

g = 0

b = 0

turtle.colormode(255)


turtle.begin_fill()

for i in range(255):

    r = (r+3) % 25

    g = (g+5) % 255

    b = (b+7) % 255

    turtle.pencolor(r,g,b)

    turtle.forward(200)

    turtle.left(360/m*n)

turtle.end_fill()


turtle.done()

pythonのトレーニングをしたと思うことにする。コマンドラインでm,nは入力しているが,簡単な有理数でないものを入れたい場合は,m=3.1415926 n=1 を代入すればよい。

図:pythonのΠ角形タートルグラフィックスの途中

2023年1月16日月曜日

Π角形

YouTubeでm/n角形という話題を見かけた。

これは,タートルグラフィックスの定番らしい。Processingでプログラムを作ってみることにする。久々のProcessingはロゴも変わりバージョンは4.1.1になっている。

タートルグラフィックスの見本コードがサンプルの中にあったので,それを参考にしてみる。クリックするとタートルが進むとして,k角形の形は,一辺の長さを進むタートルの進行方向(方向単位ベクトル)の方位角が,進むたびに360度/kづつ増加するものと定義する。

float turtleX;
float turtleY;
float turtleT;
float k;

void setup() {
  size(512, 512);
  turtleX = width/2;
  turtleY = height/2;
  turtleT = 0;
  k=3.1415926;
  background(255);
}

void forward(float go) {
  float newX = turtleX + go*cos(turtleT);
  float newY = turtleY + go*sin(turtleT);
  line(turtleX, turtleY, newX, newY);
  turtleX = newX;
  turtleY = newY;
}

void rotate(float rot) {
  turtleT = turtleT - radians(rot);
}

void mousePressed() {
  strokeWeight(2);
  stroke(mouseX/2,mouseY/2,random(0,255));
  forward(128);
  rotate(360/k);
}

ここではm/nとして簡単な有理数ではなくて,円周率の近似値(有効数字8桁)を与えてみた。つまり,Π角形の作図になる。マウスクリックでタートルが進むようにした上で,線の色はクリック位置と乱数の組み合わせで選ぶことにする。


図:Π角形作図の途中段階のイメージ

2022年12月20日火曜日

mecab-neologd(3)

mecab-neologd(2)からの続き

とりあえず,最新単語を含む形態素解析のmecab-ipadic-neologdが使えるようになった。次の課題は,これをpython プログラムの中から使えるようにすることだ。これは事例がたくさんころがっていたので,そのまま写経してみた。

2018年,Qiitaにsudo5in5kさんが書いた,mecab + NEologd + python3 で形態素解析という記事があるのでそのまま使える。import MeCab のCが小文字のtypoになっているところでつまづいた。その後,辞書ディレクトリを自分の環境に合わせて指定したところ,青空文庫の人間失格のテキストファイルをとってきて,無事に形態素分解した結果が出力された。

そこで,入出力ファイルをコマンドラインで指定し,任意のテキストファイルを変換できるように数行だけ修正したのが以下のコードである(ほぼsudo5in5kさんのもの)。

#!/opt/homebrew/bin/python3


# usage: aozora.py infile outfile

# infile is taken from https://www.aozora.gr.jp -> txt download


import MeCab

import re

import sys


infile = sys.argv[1]

outfile = sys.argv[2]


bindata = open(infile, 'rb').read()

textdata = bindata.decode('shift_jis')


# 青空文庫のための固有処理

textdata = re.split(r'\-{5,}', textdata)[2]

textdata = re.split(r'底本:', textdata)[0]

textdata = textdata.strip()


# 人によっては以下のパスは異なるので確認してね

mecab = MeCab.Tagger('-d /opt/homebrew/lib/mecab/dic/mecab-ipadic-neologd')

mecab.parse('')  # バグ対処

results = []

lines = textdata.split("\r\n")

for line in lines:

    r = []

    # 学習に使わない表現の削除処理

    s = line

    s = s.replace("|", "")

    s = re.sub(r'《.+?》', "", s)

    s = re.sub(r'[.+?]', '', s)

    # Mecab

    node = mecab.parseToNode(s)

    while node:

        # 単語を取得

        if node.feature.split(",")[6] == '*':

            word = node.surface

        else:

            word = node.feature.split(",")[6]


        # 品詞を取得

        part = node.feature.split(",")[0]


        if part in ["名詞", "形容詞", "動詞", "記号"]:

            r.append(word)

        node = node.next

    rl = (" ".join(r)).strip()

    results.append(rl)


# write to a file

with open(outfile, 'w', encoding='utf-8') as wf:

    wf.write("\n".join(results))