Python in the box – Leave record about python coding

11月 17, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST3

Each result returned by the distributed processing is not made into a PN-sized matrix, but is placed as is. This lightens the processing.

import random
import math
import time
import gc
import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

core = 6

#total particle num
PN = 10000
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 1000
xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))


def wrapper(args):
  return calc_XY_start_end(*args)

def find_pair():
  global core
  global PN
  global xyzF
  global xyzP

  worklist = []
  thread   = 0

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  Thead_num = 0
  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])
      Thead_num = Thead_num + 1
      #----------------------------------------
      #local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------


  #start thread. results is callback in array.
  p = Pool(core)
  callback = p.map(wrapper, worklist)
  p.close()


  t1 = 0
  for Th in range(0,Thead_num):

    xn_l = callback[Th][0]
    yn_l = callback[Th][1]
    xxst = callback[Th][2]
    xxed = callback[Th][3]
    yyst = callback[Th][4]
    yyed = callback[Th][5]
    local_force = callback[Th][6]

    time_sta = time.time()

    sb = local_force.shape
    ten   = -1 * local_force.transpose(1,0,2)
    tensb = ten.shape

    sum_f[yyst:yyst+sb[0], xxst:xxst+sb[1]] += local_force

    if xn_l != yn_l:
      sum_f[xxst:xxst+tensb[0], yyst:yyst+tensb[1]] += ten

    time_end = time.time()
    t1 =t1 + (time_end  - time_sta)


  print ("=====t1=",t1)
  print(np.sum(sum_f))
  #mmm = sum_f.sum(axis=1) 
  #xyzF = xyzF + mmm 


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
#  while pnum < PN:
#    total_F = total_F + xyzF[pnum][XX]
#    total_F = total_F + xyzF[pnum][YY]
#    total_F = total_F + xyzF[pnum][ZZ]
#    pnum += 1
#
#  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

import random

import math

import time

import gc

import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

core = 6

#total particle num

PN = 10000

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

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

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

divide_N = 1000

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def wrapper(args):

return calc_XY_start_end(*args)

def find_pair():

global core

global PN

global xyzF

global xyzP

worklist = []

thread = 0

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

Thead_num = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])

Thead_num = Thead_num + 1

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

#local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

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

#start thread. results is callback in array.

p = Pool(core)

callback = p.map(wrapper, worklist)

p.close()

t1 = 0

for Th in range(0,Thead_num):

xn_l = callback[Th][0]

yn_l = callback[Th][1]

xxst = callback[Th][2]

xxed = callback[Th][3]

yyst = callback[Th][4]

yyed = callback[Th][5]

local_force = callback[Th][6]

time_sta = time.time()

sb = local_force.shape

ten = -1 * local_force.transpose(1,0,2)

tensb = ten.shape

sum_f[yyst:yyst+sb[0], xxst:xxst+sb[1]] += local_force

if xn_l != yn_l:

sum_f[xxst:xxst+tensb[0], yyst:yyst+tensb[1]] += ten

time_end = time.time()

t1 =t1 + (time_end - time_sta)

print ("=====t1=",t1)

print(np.sum(sum_f))

#mmm = sum_f.sum(axis=1)

#xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

# while pnum < PN:

# total_F = total_F + xyzF[pnum][XX]

# total_F = total_F + xyzF[pnum][YY]

# total_F = total_F + xyzF[pnum][ZZ]

# pnum += 1

# print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

11月 16, 202211月 16, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST2

The split process reduced the matrix size as much as possible, but it still eats up too much memory.
Multiprocessing consumes a lot of memory.

import random
import math
import time
import gc
import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 1200
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 50
core = 2
xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))

def wrapper(args):
  return calc_XY_start_end(*args)

def find_pair():
  global PN
  global xyzF
  global xyzP

  worklist = []
  thread   = 0

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  Thead_num = 0
  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])
      Thead_num = Thead_num + 1
      #----------------------------------------
      #local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------

  p = Pool(core)

  #start thread. results is callback in array.
  callback = p.map(wrapper, worklist)
  p.close()


  for Th in range(0,Thead_num):

    xn_l = callback[Th][0]
    yn_l = callback[Th][1]
    xxst = callback[Th][2]
    xxed = callback[Th][3]
    yyst = callback[Th][4]
    yyed = callback[Th][5]
    local_force = callback[Th][6]

    if xxst != 0:
      r2 = np.zeros((xxst,divide_N,3))
      local_force        = np.insert(local_force     ,0,r2,axis=1)
    if xxed+1 != PN:
      r2 = np.zeros((PN - xxed - 1,divide_N,3))
      local_force      = np.insert(local_force     ,xxed + 1,r2,axis=1)

    if yyst != 0:
      r2 = np.zeros((yyst,PN,3))
      local_force      = np.insert(local_force     ,0,r2,axis=0)
    if yyed+1 != PN:
      r2 = np.zeros((PN - yyed - 1,PN,3))
      local_force      = np.insert(local_force     ,yyed + 1,r2,axis=0)

    del r2

    if xn_l != yn_l:
      sum_f = sum_f + local_force
      sum_f = sum_f - local_force.transpose(1,0,2)
    else:
      sum_f = sum_f + local_force

    del local_force
    gc.collect()

  #print (sum_f)

  print("=====================================")
  mmm = sum_f.sum(axis=1)
  xyzF = xyzF + mmm


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

import random

import math

import time

import gc

import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 1200

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

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

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

divide_N = 50

core = 2

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def wrapper(args):

return calc_XY_start_end(*args)

def find_pair():

global PN

global xyzF

global xyzP

worklist = []

thread = 0

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

Thead_num = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])

Thead_num = Thead_num + 1

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

#local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

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

p = Pool(core)

#start thread. results is callback in array.

callback = p.map(wrapper, worklist)

p.close()

for Th in range(0,Thead_num):

xn_l = callback[Th][0]

yn_l = callback[Th][1]

xxst = callback[Th][2]

xxed = callback[Th][3]

yyst = callback[Th][4]

yyed = callback[Th][5]

local_force = callback[Th][6]

if xxst != 0:

r2 = np.zeros((xxst,divide_N,3))

local_force = np.insert(local_force ,0,r2,axis=1)

if xxed+1 != PN:

r2 = np.zeros((PN - xxed - 1,divide_N,3))

local_force = np.insert(local_force ,xxed + 1,r2,axis=1)

if yyst != 0:

r2 = np.zeros((yyst,PN,3))

local_force = np.insert(local_force ,0,r2,axis=0)

if yyed+1 != PN:

r2 = np.zeros((PN - yyed - 1,PN,3))

local_force = np.insert(local_force ,yyed + 1,r2,axis=0)

del r2

if xn_l != yn_l:

sum_f = sum_f + local_force

sum_f = sum_f - local_force.transpose(1,0,2)

else:

sum_f = sum_f + local_force

del local_force

gc.collect()

#print (sum_f)

print("=====================================")

mmm = sum_f.sum(axis=1)

xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

11月 16, 202211月 16, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST1

The small distributed matrices also eat too much memory because they are only PNxPN in size.
Too slow.

###########################
# multiprocessing test ####
###########################
import random
import math
import time
import scipy.special as scm
from multiprocessing import Pool

random.seed(1)

PX = 0;PY = 1;PZ = 2;
VX = 3;VY = 4;VZ = 5;
FX = 6;FY = 7;FZ = 8;

#number of particles in a line
line_num = 15

#total particle num
PN = line_num * line_num * line_num

#ready to 9 parameters for particle 
#(PX, PY, PZ, VX, VY, VZ, FX, FY, FZ)
xyz = [[0 for i in range(9)] for j in range(PN)]

#Number of combinations of coulomb force calculation
combinum = int(scm.comb(PN, 2))

#thread number(local thread num)
core = 4

def find_pair_sub(prep,pend,thread):
  global xyz

  #local results array
  xyzF = [[0 for i in range(3)] for j in range(PN)]
  fx = 0; fy = 1; fz = 2

  for i in range(prep,pend):
    for j in range(i + 1, PN):
      dx = xyz[i][PX] - xyz[j][PX]
      dy = xyz[i][PY] - xyz[j][PY]
      dz = xyz[i][PZ] - xyz[j][PZ]
      r  = math.sqrt(dx*dx + dy*dy + dz*dz)

      xyzF[i][fx] = xyzF[i][fx] + dx/(r*r*r)
      xyzF[i][fy] = xyzF[i][fy] + dy/(r*r*r)
      xyzF[i][fz] = xyzF[i][fz] + dz/(r*r*r)
      xyzF[j][fx] = xyzF[j][fx] - dx/(r*r*r)
      xyzF[j][fy] = xyzF[j][fy] - dy/(r*r*r)
      xyzF[j][fz] = xyzF[j][fz] - dz/(r*r*r)

  return xyzF

def wrapper(args):
  return find_pair_sub(*args)


def find_pair():
  global PN
  global combinum

  pw = combinum // core
  pl = combinum % core

  localt = 0
  thread = 0
  pre = 0
  #each thread work list
  worklist = []
  ppp = pw

  for i in range(PN) :
    if core == 1:
      worklist.append([pre,PN,thread])
      break

    localt = localt + (PN - i - 1)
    if localt >= ppp:
      worklist.append([pre,i,thread])
      ppp += pw
      thread += 1
      pre = i

  if i != pre:
    prep = worklist[thread-1][0]
    worklist[thread-1] = [prep,PN,thread-1]

  #make thread core num
  p = Pool(core)

  #start thread. results is callback in array.
  callback = p.map(wrapper, worklist)
  p.close()

  #summation each thread results
  for j in range(core):
    for i in range(PN):
      xyz[i][FX] += callback[j][i][0]
      xyz[i][FY] += callback[j][i][1]
      xyz[i][FZ] += callback[j][i][2]

def init_lattice():
  global xyz

  pnum = 0
  while pnum < PN:
    xyz[pnum][PX] = random.uniform(-1,1)
    xyz[pnum][PY] = random.uniform(-1,1)
    xyz[pnum][PZ] = random.uniform(-1,1)
    xyz[pnum][FX] = random.uniform(-1,1)
    xyz[pnum][FY] = random.uniform(-1,1)
    xyz[pnum][FZ] = random.uniform(-1,1)
    pnum += 1

if __name__ == "__main__":
  init_lattice()
  find_pair()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

###########################

# multiprocessing test ####

###########################

import random

import math

import time

import scipy.special as scm

from multiprocessing import Pool

random.seed(1)

PX = 0;PY = 1;PZ = 2;

VX = 3;VY = 4;VZ = 5;

FX = 6;FY = 7;FZ = 8;

#number of particles in a line

line_num = 15

#total particle num

PN = line_num * line_num * line_num

#ready to 9 parameters for particle

#(PX, PY, PZ, VX, VY, VZ, FX, FY, FZ)

xyz = [[0 for i in range(9)] for j in range(PN)]

#Number of combinations of coulomb force calculation

combinum = int(scm.comb(PN, 2))

#thread number(local thread num)

core = 4

def find_pair_sub(prep,pend,thread):

global xyz

#local results array

xyzF = [[0 for i in range(3)] for j in range(PN)]

fx = 0; fy = 1; fz = 2

for i in range(prep,pend):

for j in range(i + 1, PN):

dx = xyz[i][PX] - xyz[j][PX]

dy = xyz[i][PY] - xyz[j][PY]

dz = xyz[i][PZ] - xyz[j][PZ]

r = math.sqrt(dx*dx + dy*dy + dz*dz)

xyzF[i][fx] = xyzF[i][fx] + dx/(r*r*r)

xyzF[i][fy] = xyzF[i][fy] + dy/(r*r*r)

xyzF[i][fz] = xyzF[i][fz] + dz/(r*r*r)

xyzF[j][fx] = xyzF[j][fx] - dx/(r*r*r)

xyzF[j][fy] = xyzF[j][fy] - dy/(r*r*r)

xyzF[j][fz] = xyzF[j][fz] - dz/(r*r*r)

return xyzF

def wrapper(args):

return find_pair_sub(*args)

def find_pair():

global PN

global combinum

pw = combinum // core

pl = combinum % core

localt = 0

thread = 0

pre = 0

#each thread work list

worklist = []

ppp = pw

for i in range(PN) :

if core == 1:

worklist.append([pre,PN,thread])

break

localt = localt + (PN - i - 1)

if localt >= ppp:

worklist.append([pre,i,thread])

ppp += pw

thread += 1

pre = i

if i != pre:

prep = worklist[thread-1][0]

worklist[thread-1] = [prep,PN,thread-1]

#make thread core num

p = Pool(core)

#start thread. results is callback in array.

callback = p.map(wrapper, worklist)

p.close()

#summation each thread results

for j in range(core):

for i in range(PN):

xyz[i][FX] += callback[j][i][0]

xyz[i][FY] += callback[j][i][1]

xyz[i][FZ] += callback[j][i][2]

def init_lattice():

global xyz

pnum = 0

while pnum < PN:

xyz[pnum][PX] = random.uniform(-1,1)

xyz[pnum][PY] = random.uniform(-1,1)

xyz[pnum][PZ] = random.uniform(-1,1)

xyz[pnum][FX] = random.uniform(-1,1)

xyz[pnum][FY] = random.uniform(-1,1)

xyz[pnum][FZ] = random.uniform(-1,1)

pnum += 1

if __name__ == "__main__":

init_lattice()

find_pair()

11月 15, 2022

The process of brute force calculation by Coulomb force

The process of brute force calculation by Coulomb force is calculated as a matrix using numpy. In the calculation, each element is calculated in a small matrix and summed up at the end to enable parallel processing later.

import random
import math
import time
import numpy as np

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 100
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 3

xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))

def find_pair():
  global PN
  global xyzF
  global xyzP

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      #----------------------------------------
      local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------

      if Xn != Yn:
        sum_f = sum_f + local_force
        sum_f = sum_f - local_force.transpose(1,0,2)
      else:
        sum_f = sum_f + local_force

  print (sum_f)

  print("=====================================")
  mmm = sum_f.sum(axis=1)
  xyzF = xyzF + mmm


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  if X_st != 0:
    r2 = np.zeros((X_st,divide_N,3))
    tmp_f        = np.insert(tmp_f     ,0,r2,axis=1)
  if X_ed+1 != PN:
    r2 = np.zeros((PN - X_ed - 1,divide_N,3))
    tmp_f      = np.insert(tmp_f     ,X_ed + 1,r2,axis=1)

  if Y_st != 0:
    r2 = np.zeros((Y_st,PN,3))
    tmp_f      = np.insert(tmp_f     ,0,r2,axis=0)
  if Y_ed+1 != PN:
    r2 = np.zeros((PN - Y_ed - 1,PN,3))
    tmp_f      = np.insert(tmp_f     ,Y_ed + 1,r2,axis=0)

  #print (tmp_f)
  #print("===================================================")
  return tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

import random

import math

import time

import numpy as np

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 100

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

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

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

divide_N = 3

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def find_pair():

global PN

global xyzF

global xyzP

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

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

local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

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

if Xn != Yn:

sum_f = sum_f + local_force

sum_f = sum_f - local_force.transpose(1,0,2)

else:

sum_f = sum_f + local_force

print (sum_f)

print("=====================================")

mmm = sum_f.sum(axis=1)

xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

if X_st != 0:

r2 = np.zeros((X_st,divide_N,3))

tmp_f = np.insert(tmp_f ,0,r2,axis=1)

if X_ed+1 != PN:

r2 = np.zeros((PN - X_ed - 1,divide_N,3))

tmp_f = np.insert(tmp_f ,X_ed + 1,r2,axis=1)

if Y_st != 0:

r2 = np.zeros((Y_st,PN,3))

tmp_f = np.insert(tmp_f ,0,r2,axis=0)

if Y_ed+1 != PN:

r2 = np.zeros((PN - Y_ed - 1,PN,3))

tmp_f = np.insert(tmp_f ,Y_ed + 1,r2,axis=0)

#print (tmp_f)

#print("===================================================")

return tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

9月 29, 2022

Calculating long-range particle interactions with numpy (handmade mesh grid)

Too slow

import random
import math
import time
import numpy as np

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 300

xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))


def find_pair():
  global PN
  global xyzF
  global xyzP

  x1 = np.empty(0)
  y1 = np.empty(0)

  for mm in range(1,PN):
    for nn in range(mm,PN):
      x1 = np.append(x1,[nn],axis=0)
      y1 = np.append(y1,[mm-1],axis=0)

  x1 = x1.astype(np.int32)
  x1 = x1.reshape(x1.shape[0],1)
  y1 = y1.astype(np.int32)
  y1 = y1.reshape(y1.shape[0],1)

  dd = np.linalg.norm(xyzP[x1]-xyzP[y1], axis=2)

  distances = np.zeros((PN,PN))
  ly = 0
  lx = ly + 1
  for kk in range(0,dd.shape[0]):
    distances[lx][ly] = dd[kk]
    distances[ly][lx] = dd[kk]
    lx = lx + 1
    if lx == PN:
      ly = ly + 1
      lx = ly + 1

  tmp_index = np.arange(xyzP.shape[0])
  xx, yy = np.meshgrid(tmp_index, tmp_index)

  dxyz = xyzP[xx] - xyzP[yy]

  distances = distances[:,:,np.newaxis]
  tmp_force = dxyz/distances
  tmp_force = np.nan_to_num(tmp_force,0)

  mmm = tmp_force.sum(axis=1)
  xyzF = xyzF + mmm

def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1

def results_sum():
  #global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

import random

import math

import time

import numpy as np

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 300

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def find_pair():

global PN

global xyzF

global xyzP

x1 = np.empty(0)

y1 = np.empty(0)

for mm in range(1,PN):

for nn in range(mm,PN):

x1 = np.append(x1,[nn],axis=0)

y1 = np.append(y1,[mm-1],axis=0)

x1 = x1.astype(np.int32)

x1 = x1.reshape(x1.shape[0],1)

y1 = y1.astype(np.int32)

y1 = y1.reshape(y1.shape[0],1)

dd = np.linalg.norm(xyzP[x1]-xyzP[y1], axis=2)

distances = np.zeros((PN,PN))

ly = 0

lx = ly + 1

for kk in range(0,dd.shape[0]):

distances[lx][ly] = dd[kk]

distances[ly][lx] = dd[kk]

lx = lx + 1

if lx == PN:

ly = ly + 1

lx = ly + 1

tmp_index = np.arange(xyzP.shape[0])

xx, yy = np.meshgrid(tmp_index, tmp_index)

dxyz = xyzP[xx] - xyzP[yy]

distances = distances[:,:,np.newaxis]

tmp_force = dxyz/distances

tmp_force = np.nan_to_num(tmp_force,0)

mmm = tmp_force.sum(axis=1)

xyzF = xyzF + mmm

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

#global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

9月 28, 2022

Calculating long-range particle interactions with numpy

import random
import math
import time
import numpy as np

#import scipy.misc as scm
from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 1000

xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))

def find_pair():
  global PN
  global xyzF
  global xyzP

  tmp_index = np.arange(xyzP.shape[0])

  xx, yy = np.meshgrid(tmp_index, tmp_index)
  distances = np.linalg.norm(xyzP[xx]-xyzP[yy], axis=2)
  dxyz = xyzP[xx] - xyzP[yy]

  distances = distances[:,:,np.newaxis]
  tmp_force = dxyz/distances
  tmp_force = np.nan_to_num(tmp_force,0)

  mmm = tmp_force.sum(axis=1)
  xyzF = xyzF + mmm

def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  #global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

import random

import math

import time

import numpy as np

#import scipy.misc as scm

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 1000

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def find_pair():

global PN

global xyzF

global xyzP

tmp_index = np.arange(xyzP.shape[0])

xx, yy = np.meshgrid(tmp_index, tmp_index)

distances = np.linalg.norm(xyzP[xx]-xyzP[yy], axis=2)

dxyz = xyzP[xx] - xyzP[yy]

distances = distances[:,:,np.newaxis]

tmp_force = dxyz/distances

tmp_force = np.nan_to_num(tmp_force,0)

mmm = tmp_force.sum(axis=1)

xyzF = xyzF + mmm

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

#global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

9月 25, 20229月 25, 2022

scipy.misc を scipy.special に置き換え

AttributeError: module ‘scipy.misc’ has no attribute ‘comb’

というエラーが出るようになった。

combを使うに当たっては、scipy.special から呼び出さなくてはいけない様子。

https://python-forum.io/thread-21470.html
https://docs.scipy.org/doc/scipy-1.2.1/reference/misc.html

import scipy.misc
↓
import scipy.special

の書き換えで解決。

1月 14, 20221月 14, 2022

乱数を使って円周率もどきを求める2

先日試した乱数を使った円周率を求める方法を
numpyを使った方法に置き換えた。

import math
import time
import numpy as np


def calc_pi ( NN ):
  arr1 = np.random.rand(NN)
  arr2 = np.random.rand(NN)
  arr3 = np.sqrt(arr1 * arr1 + arr2 * arr2)
  arr3[arr3 > 1] = 0   #1より大きい要素は0に
  arr3[arr3 > 0] = 1   #0より大きい要素は1に
  in_n = np.sum(arr3)  #円内に入ったやつの合計

  return ( in_n / NN ) * 4


if __name__ == "__main__":

  for j in range (1,9):
    NN = 10 ** j
    start_time = time.perf_counter()
    pi = calc_pi(NN)
    execution_time = time.perf_counter() - start_time
    delt = math.fabs(math.pi - pi)
    print (str(NN).rjust(10),"  ",f'{pi:11.010f}',"  ",f'{delt:7.06f}',"  ",execution_time)

import math

import time

import numpy as np

def calc_pi ( NN ):

arr1 = np.random.rand(NN)

arr2 = np.random.rand(NN)

arr3 = np.sqrt(arr1 * arr1 + arr2 * arr2)

arr3[arr3 > 1] = 0 #1より大きい要素は0に

arr3[arr3 > 0] = 1 #0より大きい要素は1に

in_n = np.sum(arr3) #円内に入ったやつの合計

return ( in_n / NN ) * 4

if __name__ == "__main__":

for j in range (1,9):

NN = 10 ** j

start_time = time.perf_counter()

pi = calc_pi(NN)

execution_time = time.perf_counter() - start_time

delt = math.fabs(math.pi - pi)

print (str(NN).rjust(10)," ",f'{pi:11.010f}'," ",f'{delt:7.06f}'," ",execution_time)

実行結果

試行回数結果　　　　　　誤差　　　　　処理時間[sec]

10 2.8000000000 0.341593 8.769400028540986e-05

100 3.1200000000 0.021593 4.6955000016168924e-05

1000 3.1960000000 0.054407 5.6955000218295027e-05

10000 3.1420000000 0.000407 0.0003660469997157634

100000 3.1376400000 0.003953 0.0036325510000096983

1000000 3.1415240000 0.000069 0.03182136600025842

10000000 3.1416492000 0.000057 0.2625446400002147

100000000 3.1416174400 0.000025 2.535640128000068

numpy使わずにwhileで計算していたまえの処理に較べて15倍位早い結果。全て配列に収めてから処理するとさすがに早い。ただしメモリもバカ喰い。

1月 13, 20221月 14, 2022

乱数を使って円周率もどきを求める

乱数を使った円周率もどきを求めてみる

下記のような正方形と内接する円を考える。正方形の1辺は２としておく。この正方形の中にランダムに点を打った際、円の中にも入る確率は円の面積 / 正方形の面積。つまり、確率 = π / 4 となる。

π = 正方形内にランダムに点を打った際に円の中に入る確率は × 4

ランダムに打った点が円の中に入ったかどうか？はｘ座標とｙ座標にそれぞれ0～1の乱数を与えてやり、sqrt(x^2 + y^2) が1以下であれば円の中といえる。

この処理を実装したものが下記

import random
import math
import time

def calc_rad (x, y):
  return math.sqrt(x*x +  y*y)


def calc_pi ( NN ):
  num  = 0
  in_n = 0
  while num < NN:
    x = random.random()
    y = random.random()
    if calc_rad (x, y) <= 1:
      in_n = in_n  + 1
    num = num + 1

  return ( in_n / num ) * 4


if __name__ == "__main__":
  for j in range (1,9):
    NN = 10 ** j
    start_time = time.perf_counter()
    pi = calc_pi(NN)
    execution_time = time.perf_counter() - start_time
    delt = math.fabs(math.pi - pi)
    print (str(NN).rjust(10),"  ",f'{pi:11.010f}',"  ",f'{delt:7.06f}',"  ",execution_time)

import random

import math

import time

def calc_rad (x, y):

return math.sqrt(x*x + y*y)

def calc_pi ( NN ):

num = 0

in_n = 0

while num < NN:

x = random.random()

y = random.random()

if calc_rad (x, y) <= 1:

in_n = in_n + 1

num = num + 1

return ( in_n / num ) * 4

if __name__ == "__main__":

for j in range (1,9):

NN = 10 ** j

start_time = time.perf_counter()

pi = calc_pi(NN)

execution_time = time.perf_counter() - start_time

delt = math.fabs(math.pi - pi)

print (str(NN).rjust(10)," ",f'{pi:11.010f}'," ",f'{delt:7.06f}'," ",execution_time)

実行結果

試行回数結果　　　　　　誤差　　　　　処理時間[sec]

10 3.6000000000 0.458407 1.1099999937869143e-05

100 2.9600000000 0.181593 4.040199928567745e-05

1000 3.0800000000 0.061593 0.0004007130000900361

10000 3.1616000000 0.020007 0.004138636999414302

100000 3.1487600000 0.007167 0.042758481000419124

1000000 3.1437800000 0.002187 0.4172533930004647

10000000 3.1411992000 0.000393 4.038885087000381

100000000 3.1415664400 0.000026 41.15365210999971

それなりに近い値かな？

1月 10, 20221月 10, 2022

centos7にてpycudaをインストール(pip3)

pycudaを下記コマンドからインストールを試したところ、、、

1 2	pip3 install pycuda

エラーでストップ。
pycudaのコンパイル時にcuda.hがないとのメッセージで死んでいる。

このエラーについてググるとコンパイル時のPATH/
CPATH/LIBRARY_PATH等にcudaのpathを追加で
対応可能と出てくるが自分の環境では改善せず。。

自分のところでは、gccのバージョンを上げることでエラーが消えて
インストール成功した。gccの対応バージョンがあるっぽい。

エラー時：gcc version 4.8.5 20150623
成功時　：gcc version 7.3.1 20180303