Using pyplot@Plots on Julia@spack

 Using pyplot@Plots on Julia installed via spack is a bit complicated.

We first need to setup Python environment
(actually, you may want to use conda provided by Julia, so, you may do
Julia> ENV["PYTHON"]=""
)

as described in the following link;
https://github.com/JuliaPy/PyCall.jl

and then,

Julia> using Plots
Julia> pyplot()

Normal distribution of dice rolls: Python animation

 Since my wife started learning statistics, and I wanted to give her some intuitive understanding of the normal distribution, standard deviation, 95% confidence region etc, I made a small python code which generates an animation of the histogram of dice rolls.

This could help others time ;)

(Distribute under MIT license)

import random

import matplotlib.pyplot as plt

import numpy as np

from scipy.stats import norm

def main():

numTrial = 1000000

numDice = 20

sumResults = []

bin=np.arange(numDice-0.5, numDice*6+1.5, 1)

print(bin)

fig = plt.figure()

iFig = 1

for iTrial in range(numTrial):

sumDice = 0

for iDice in range(numDice):

sumDice += random.randint(1,6)

sumResults.append(sumDice)

if iTrial % 100 == 0:

ax = fig.add_subplot()

ax.set_title(str(numDice)+" dices case")

ax.set_xlabel("Sum of rolls")

ax.set_ylabel("Frequency")

n, bins, patches = plt.hist(sumResults,bins=bin)

loc,scale = norm.fit(sumResults)

plt.plot(bin,sum(n)*norm.pdf(bin,loc,scale),color="y")

txt = "STD:"+str(round(scale,4))

plt.text(5.0,0.01,txt,size=30)

p25 =norm.ppf(0.025 ,loc,scale)

p975=norm.ppf(0.975,loc,scale)

plt.axvline(x=p25, ymin=0.1, ymax=0.9,color="r")

plt.axvline(x=p975, ymin=0.1, ymax=0.9,color="r")

plt.savefig(str('{:05d}'.format(iFig))+".png")

iFig += 1

plt.pause(.001)

fig.clear()

if __name__ == "__main__":

main()

Colour scheme for scientific computing visualization

 Looks nice, and the author provides Python codes (and some other languages).

I also have a kind of colour blindness (red and green). It would be nice to have this.

https://ai.googleblog.com/2019/08/turbo-improved-rainbow-colormap-for.html

Data transfer from Oracle To PostgreSQL via Pandas

 I needed to transfer data from Oracle to PostgreSQL, which are on different machines.
This time, I used Pandas and Pickle on Python 2.7.

Reading from Oracle:

import pickle

import cx_Oracle

import pandas as pd

def main():

conn = cx_Oracle.connect("ID/PW@DB")

tables = ["a","b"]

data = {}

for iTab in tables:

query = "select * from " + iTab

curDF = pd.read_sql(query, con=conn)

data[iTab] = curDF

conn.close()

with open("data.dat","wb") as f:

pickle.dump(data,f)

Writing to PostgreSQL:
(Note1, table name, schema and column names should be in lower case. Otherwise, you will find or need strange double quotation marks.
Note2, you should give chunksize. Otherwise, the python instance will eat up all the memory.
Note3, you need to set up schemas, but tables will be automatically created)

import pickle

import pandas as pd

from sqlalchemy import create_engine

def main():

with open("data.dat","rb") as f:

arrData = pickle.load(f)

engine = create_engine('postgresql://ID:PW@localhost:5432/server')

for iTab in arrData:

print(iTab)

arrData[iTab].to_sql(iTab.lower(), engine,schema="schema".lower(),chunksize=100)

Jug, Parallel framework in Python

 Minor tips for problems I was into in Jug, a parallel computing framework in Python.
I may be wrong, and if you know I misunderstand, please leave comments:

1. Tasks can be created only outside functions.
   Probably, this is common to create the main function with, if __name__ = __main__: main()
   I did this, and called functions decorated with @TaskGenerator. However, jug could not launch tasks that are created in the main function. Instead, when I remove the main function, and write lines on the ground level (no indent level), the same code worked.
2. Multi-node parallelism can be achieved by just running "jug execute XXX.py" if the directory which stores the python code is shared by NFS. I could not find any clear comment on this, and I wondered how I can achieve multi-node parallelism.
3. An atomic operation for File I/O can be achieved by creating a function decorated with @TaskGenerator. According to the summary of Jug, such a function is executed once. So, we can use it for an atomic operation.

Otherwise, I think Jug is excellent, and easy to use!

Still a question:

I do not know what jugfile can do.
On page17 of the documentation(https://jug.readthedocs.io/_/downloads/en/latest/pdf/), a post-process script import jugfile. I think the tasks in the main script make outputs to the jugfile, but do not know how to do.

Upper and lower limit in matplotlib contourf

When the data range in matplotlib contourf does not correspond to the range of colorbar, matplotlib automatically adjust the range, even when set_clim is set manually.

(In other words, I needed a fixed colorbar regardless of the data)

The solution is set "levels" as follows;

interval = np.linspace(2.5,4.5, num=32, endpoint=True)

img=plt.contourf(LON, LAT, DATA,

vmin=2.5,vmax=4.5,

cmap = "jet",

levels=interval,

transform=ccrs.PlateCarree())

Pandas on Python 2.6

We are still using Python2.6 on RH6, and packages are getting out-dated.

Pandas is one of those, and need the following command to be installed.

% easy_install 'pandas==0.17.1'

This seems the latest version for Python 2.6

And also numexpr is necessary for the "query" but the latest version numexpr cannot be installed.

% easy_install 'pandas==2.4'

at least works.

Convert datetime to epoch in a old Python (2.6)

 Since I always forget how to convert datetime to epoch, I just want to write a memo.

In Python
>>> import datetime
>>> d = datetime.datetime(2019,1,1,0,0,0)
>>> print(d)

2019-01-01 00:00:00

>>> epoch = d.strftime("%s")
>>> print(epoch)

1546300800

On Bash

% date -d @1546300800

Tue Jan  1 00:00:00 GMT 2019

Leetcode: Unique Paths

 To make some practice on DP, "Unique Paths"@Leetcode is solved.

Probably a typical one, and super easy.
https://leetcode.com/problems/unique-paths

import numpy as np

class Solution:

def uniquePaths(self, m: int, n: int) -> int:

dp = np.zeros((m,n),np.int32)

dp[0,:]=1

dp[:,0]=1

for j in range(1,m):

for i in range(1,n):

dp[j,i] = dp[j-1,i] + dp[j,i-1]

return dp[m-1][n-1]

Sudoku solver @LeetCode

 Not elegant at all, but readable sudoku solver.

https://leetcode.com/problems/sudoku-solver/

class Solution(object):

    

    def checkNew(self, i,j,k,board):

        

        for jj in range(9):

            if board[i][jj] == str(k):

                return False

        for ii in range(9):

            if board[ii][j] == str(k):

                return False

        

        iB = i//3

        jB = j//3

        b = []

        for ii in range(3):

            for jj in range(3):

                if board[3*iB+ii][3*jB+jj] == str(k):

                    return False

        

        

        return True

    

    def backtrack(self, board):

# Find a blank element. The search direction does not a matter,

# but row-oriented here.

        for i in range(9):

            for j in range(9):

                if board[i][j] == ".":

# If blank, try numbers from 1 to 9

                    for k in range(9):

                        if self.checkNew(i,j,k+1,board): # Check if k is OK            

                            board[i][j]=str(k+1)

                            if self.backtrack(board): # Try next blank until OK

                                return True

                            else:

                                board[i][j] = "." # if inconsistent, step back

                    else:

                        return False

        return True

    

    

    def solveSudoku(self, board):

        """

        :type board: List[List[str]]

        :rtype: None Do not return anything, modify board in-place instead.

        """

        self.backtrack(board)

Unique Paths III @LeetCode

Just to wake up; 

980. Unique Paths III
https://leetcode.com/problems/unique-paths-iii/

class Solution:

    

    def checkIfReachable(self,curPoint,walkable, grid):

        (i,j) = curPoint

        #print(curPoint)

        if i < 0 or j < 0 or i >= len(grid) or j >= len(grid[0]):

            return 0

        

        if grid[i][j]==2:

            if len(walkable)==0:

                return 1

            else:

                return 0

        

        if curPoint not in walkable:

            return 0

        

        walkable.remove(curPoint)

        

        num1 = self.checkIfReachable((i+1,j  ),walkable.copy(),grid)

        num2 = self.checkIfReachable((i-1,j  ),walkable.copy(),grid)

        num3 = self.checkIfReachable((i  ,j+1),walkable.copy(),grid)

        num4 = self.checkIfReachable((i  ,j-1),walkable.copy(),grid)

        

        return num1+num2+num3+num4

    

    def uniquePathsIII(self, grid: List[List[int]]) -> int:

        

        walkable = []

        

        for i in range(len(grid)):

            for j in range(len(grid[0])):

                if grid[i][j] == 0:

                    walkable.append((i,j))

                elif grid[i][j] == 1:

                    walkable.append((i,j))

                    start = (i,j)

        

        numPath = self.checkIfReachable(start,walkable,grid)

        return numPath

TDPC C

 というわけで、TDPC　C。

考え方は簡単で、
(1) 自分の対戦相手が今まで生き残ってくる確率
(2) 自分が勝てる確率
(3) 自分が今まで生き残って来れた確率

を掛け算して足しあわせると現在自分が勝ち上がれるかどうかの確率になる
（かどうか、いまいち直感的でないんだけど、答えはそれであってる）。

どっちかというと対戦相手を見つけてくるループが面倒で、私は予めリストアップした(groupがそれ)

https://atcoder.jp/contests/tdpc/submissions/18020128

import numpy as np

def calcProb(P,Q,R):

   RP = R[P]

   RQ = R[Q]

   prob = 1.0 / (1.0+10.0**((RQ-RP)/400.0))

   return prob

def calcConsequtiveProb(iP,iG,R,dp,iBat):

   prob = 0.0

   for i in iG:

      prob += calcProb(iP,i,R)*dp[iBat-1,i]

   return prob

N = int(input())

numPerson = 2**N

R = []

for i in range(numPerson):

   R.append(float(input()))

group = []

group.append([ [i] for i in range(numPerson)])

for i in range(1,N+1):

   locGroup = []

   prevGroup = group[i-1]

   for j in range(len(prevGroup)//2):

      lg1 = prevGroup[2*j  ]

      lg2 = prevGroup[2*j+1]

      locGroup.append(lg1+lg2)

   group.append(locGroup)

dp = np.zeros((N+1,numPerson),np.float)

dp[0,:] = 1.0

for iBat in range(1,N+1):

   locGrp = group[iBat-1]

   for iGrp in range(len(locGrp)//2):

      lg1 = locGrp[iGrp*2  ]

      lg2 = locGrp[iGrp*2+1]

      for iPerson in lg1:

         prob = calcConsequtiveProb(iPerson,lg2,R,dp,iBat)

         dp[iBat,iPerson] = dp[iBat-1,iPerson]*prob

      for iPerson in lg2:

         prob = calcConsequtiveProb(iPerson,lg1,R,dp,iBat)

         dp[iBat,iPerson] = dp[iBat-1,iPerson]*prob

for i in range(numPerson):

   print(f'{dp[N,i]:.8f}')

TDPC(Typical DP Contest) A

 二次元配列にしないで解いてる人もいるみたいだけど私にはさっぱりわからないのでとりあえず二次元で。
”これまでに作成可能な点数であることが確認できれば、現在確認中のスコアを足した点数も作成可能である”ことを二次元配列で表現。

最終的には、最後の行を見れば作れる点数の個数がわかる。

https://atcoder.jp/contests/tdpc/tasks/tdpc_contest

import numpy as np

N = int(input())

line = input().split()

p = []

p.append(0)

for i in line:

   p.append(int(i))

dp = np.zeros((len(p),sum(p)+1),np.int8)

dp[0,0]=1

for j in range(1,len(p)):

   for i in range(0,sum(p)+1):

      if dp[j-1,i]:

         dp[j,i+p[j]] = 1

         dp[j,i] = 1

#print(dp)

print(sum(dp[len(p)-1,:]))

ABC129 C

DP の勉強のために下記ブログの問題リストを解いてみている。
https://qiita.com/drken/items/dc53c683d6de8aeacf5a

DPだとわかっているからできるけど、自分でいきなりはまだ思いつけなさそう。
https://atcoder.jp/contests/abc129/tasks/abc129_c

import sys

line = input().split()

N = int(line[0])

M = int(line[1])

dp = [0] * (N+1)

NA = [1] * (N+1)

for i in range(M):

   NA[int(input())] = 0

for i in range(N):

   if NA[i+1] == 0 and NA[i]==0:

      print(0)

      sys.exit(0)

dp[0] = 1

dp[1] = 1

for i in range(2,N+1):

   if NA[i] == 0:

      dp[i] = dp[i-1]

   elif NA[i-1] == 0:

      dp[i] = dp[i-1]

   else:

      dp[i] = dp[i-1]*NA[i-1] + dp[i-2]*NA[i-2]

print(dp[N]%1000000007)

 

matplotlib pcolormeshで、描画領域外でも白く塗りつぶしてしまう問題

Yin-Yang  格子の計算結果をmatplotlibで描画しているのだけれど、
Python3ではデータのない領域でも塗りつぶしてしまうようになったようで、
Yangを描画させると先に出力させていたYinの画像が塗りつぶされてしまう。
下記のURLを参考にして、データが0.0のところは透明にすることで回避できた。

領域外だからゼロじゃなくてNaNとかではないかと思うけれど、深くは気にしない。

https://teratail.com/questions/95800

ABC 158 C

S = input().split()
[A,B] = [int(i) for i in S]

Alow = A*100//8
Ahi  = (A+1)*100//8

Blow = B*100//10
Bhi  = (B+1)*100//10

#print(Alow,Ahi,Blow,Bhi)

Aran = [i for i in range(Alow,Ahi+1) if int(i*0.08) == A]
Bran = [i for i in range(Blow,Bhi+1) if int(i*0.10) == B]


resRan = sorted(list(set(Aran) & set(Bran)))
#print(resRan)
if len(resRan) !=0:
  print(resRan[0])
else:
  print(-1)

S = input().split()
[A,B] = [int(i) for i in S]

Alow = A*100//8
Ahi  = (A+1)*100//8

Blow = B*100//10
Bhi  = (B+1)*100//10

#print(Alow,Ahi,Blow,Bhi)

Aran = [i for i in range(Alow,Ahi+1) if int(i*0.08) == A]
Bran = [i for i in range(Blow,Bhi+1) if int(i*0.10) == B]


resRan = sorted(list(set(Aran) & set(Bran)))
#print(resRan)
if len(resRan) !=0:
  print(resRan[0])
else:
  print(-1)

うーん。焦ってるとこれぐらいのことも思いつかないもんだなぁ。

S = input().split()

[A,B] = [int(i) for i in S]

Alow = A*100//8

Ahi  = (A+1)*100//8

Blow = B*100//10

Bhi  = (B+1)*100//10

Aran = [i for i in range(Alow,Ahi+1) if int(i*0.08) == A]

Bran = [i for i in range(Blow,Bhi+1) if int(i*0.10) == B]

resRan = sorted(list(set(Aran) & set(Bran)))

if len(resRan) !=0:

  print(resRan[0])

else:

  print(-1)

Atcoder ABC157 E

下記コード、pypyだと通るけど、python3だとTLEする。

https://atcoder.jp/contests/abc157/submissions/10530007

from collections import defaultdict
import bisect
alphabet="abcdefghijklmnopqrstuvwxyz"
bS = int(input())
S= input()
nQ = int(input())
S = [i for i in S]
query = []
for i in range(nQ):
  query.append(input().split())
locs =defaultdict(list)
for i,char in enumerate(S):
    locs[char].append(i)
   
def countChars(iter,ss,ee):
  nChar = 0
  for i in alphabet:
    if len(locs[i]) == 0:
      pass
    else:
      sptr = bisect.bisect_right(locs[i],ss)
      eptr = bisect.bisect_left(locs[i],ee)
      if sptr == eptr:
        if sptr < 1 or sptr > len(locs[i]):
          continue
        else:
          if locs[i][sptr-1] == ss:
            nChar+=1
      else:
        nChar += 1
  return nChar
   
 
res = []
for iter,iQ in enumerate(query):
  if int(iQ[0]) == 2:
    res.append(str( countChars(iter,int(iQ[1])-1,int(iQ[2])) ))
  if int(iQ[0]) == 1:
    location = int(iQ[1])-1
    curChar  = S[location]
    if curChar != iQ[2]:
      locs[curChar].remove(location) #remove from old
      bisect.insort(locs[iQ[2]],location)
      S[location] = iQ[2]
     
print(" ".join(res))
  

pythonのソースを整形 (autopep8)

私自身はviとかEmacsでコードを書くせいで、pythonのインデントのスペースの個数が時々間違ってるらしく、eclipseを使ってる若い人に怒られたので自動整形ソフトを探す。

autopep8が簡単でうまく行った。

https://github.com/hhatto/autopep8

defaultdict

pythonで、list = [(a,1.0),(a,2.0),(b,3.0),(c,4.0)] のようなデータがあったときに、

a = [1.0, 2.0]
b = [3.0]
c = [4.0]
と整形したい。

from collections import defaultdict を使うときれいにかけて、

dict = defaultdic(list)
for i in list:
   dict[i[0]].append(i[1])

とやると、a, b, cがkeyになった辞書ができあがる。
naiveにやると二重ループになってlistが巨大なときに大変だし予めソートかけておくのも回りくどいので。

Miniconda(python2.7) + mod_wsgi on Centos6

Because the default python on Centos6 is too old, I need to use python2.7 for mapproxy.
so, I install miniconda2.
Mapproxy should work with apache, but using mod_wsgi with miniconda2 is a bit complicated.

I am describing a rough procedure:

% cd /opt/mapproxy
%  wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh

% sh Miniconda2-latest-Linux-x86_64.sh
(answer questions, “yes”, “/opt/mapproxy/miniconda2”, “no”)

% miniconda2/bin/pip install mod_wsgi
% miniconda2/bin/pip install mapproxy

% miniconda2/bin/mod_wsgi-express install-module

% vi /etc/httpd/conf/httpd.conf

LoadModule wsgi_module modules/mod_wsgi-py27.so

% chmod –R 755 /opt/mapproxy/miniconda2/

Some other configs:

WSGIScriptAlias /mapproxy /opt/mapproxy/config.py

WSGIPythonHome /opt/mapproxy/miniconda2/

WSGIApplicationGroup %{GLOBAL}

WSGISocketPrefix /var/run/wsgi

WSGIDaemonProcess mysite python-path=/opt/mapproxy/miniconda2/lib/python2.7/site-packages