흑백과 디지털

 오늘 회사에서 이것과 씨름을 했다. http://openapi.nsdi.go.kr/nsdi/eios/ServiceDetail.do?svcSe=F&svcId=F018&provOrg=NIDO 에 광역시나 도별로 GIS 건물정보를 받을 수 있다. 파일을 다운 받아서 압축을 해제하면 .shp, .dbf, .prj 파일들이 있다. 
 이 파일의 의미를 처음으로 알게되었는데,  .shp 파일은 맵의 위치, 도형정보가 들어있다. 이 파일을 잘 이용하면 실제 위도, 경도를 계산할 수 있다. .dbf 는 DB파일이다. 다운받은파일에서는 건물 주소, 건물명, 사용승인일등이 들어 있다. .shp 에 매핑되는 지점의 정보가 들어 있는 것이다. 그리고 .prj 는 .shp 파일의 맵의 위치정보가 어떤 좌표계를 사용했는지 알려주는 것 같은데 실제 이용되는 데이터인지는 모르겠다. 

 이 작업을 하면서 알게 되었는데, 우리가 흔히 위도 경로라는 것도 좌표계의 하나로 WGS84 경위도: GPS가 사용하는 좌표계를 이용해야 한다. 그리고 해당 파일의 좌표계는 "EPSG:5174" 이다. 

파이썬 코드를 작성해 보면 아래와 같다.

import shapefile
from shapely.geometry import Polygon
from dbfread import DBF
from pyproj import CRS, Transformer

# https://www.osgeo.kr/17
# EPSG:5174
EPSG5174 = """+proj=tmerc +lat_0=38 +lon_0=127.0028902777778 +k=1
+x_0=200000 +y_0=500000 +ellps=bessel +units=m +no_defs
+towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43"""


# Projection 정의
# EPSG:5174
epsg5174_crs = CRS(EPSG5174)
  

# WGS84 경위도: GPS가 사용하는 좌표계 EPSG:4326
epsg4326_crs = CRS('epsg:4326')

transformer = Transformer.from_crs(epsg5174_crs, epsg4326_crs)

fileName = './res/AL_11_D164_20230125'
sf = shapefile.Reader(shp=(fileName + ".shp"))
dbf = DBF(fileName + ".dbf" , encoding="CP949")
# records = sf.records()

index = -1
shapes = sf.shapes()
for record in dbf:
	index += 1
	if index % 10000 != 0:
		continue
	print("index : " + str(index))
	print("n_count : " + str(n_count))

	shape = shapes[index]
	polygon = Polygon(shape.points)
	center = polygon.centroid
	gis = transformer.transform(center.x, center.y)
	# gis = proj.transform(center.x, center.y)
	print("중심좌표 : " + gis)
  

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

사용한 파이썬 라이브러리는 shapefile, dbfread, pyproj 이다. shapefile은 shp 파일을 불러오는 역할을 한다. 이 라이브러리을 통해 dbf 파일도 같이 읽을 수는 있는데 속도가 느려서 dbfread 라이브러리를 사용했다. pyproj 는 shapefile 라이브러리를 통해 얻어온 좌표를 변환할 때 사용한다. 이 라이브러리를 통해 위도 경도 값을 얻을 수 있다. 
 gis 변수에 위도, 경도 tuple 이 있고, record에는 dbf 파일의 정보가 있다. 서로 같은 순서를 가지기 때문에 이렇게 하면 동일한 건물에 대해 위도,경도와 기타 정보를 같이 획득할 수 있다. 


내 설명이 많이 불친절하긴 한데, 다른 사이트들에 개념은 잘 설명되어 있다. 그런데 파이썬 코드가 없어서 이 글을 올린다.  

요즘이라고 적었지만 내 기준으로 학교에서 배우지 못한 sorting 을 정리한다. 

  Bitonic sorter(https://en.wikipedia.org/wiki/Bitonic_sorter)
 * 우선 이 sort 는 병렬 처리에 최적화된 sort 이다. 이 sort 에 대해 이해하려면 Sorting network(https://en.wikipedia.org/wiki/Sorting_network) 을 이해해야한다. Sorting Network 는 N개의 항목에 대해 비교연산을 최적화해서 sorting 하는 방법이다. 

왼쪽은 input 된 인자의 array 이고 오른쪽은 output의 array 라고 생각하면된다. sorting network 에서 위에서 위, 아래를 잇는 수직선은 두 element 끼리 비교연산을 해서 위치를 유지하거나 swap 한다고 생각하면된다. 대충 함수로 나타나면 아래와 같다. (함수에 오류가 있을 수 있다. 대충 작성했다.)

def cmp_swap(arr, idx_0, idx_1):
    if arr[idx_0] < arr[idx_1]:
        temp = arr[idx_0]
        arr[idx_1] = temp
        return arr
    return arr

def sortnetwork_desc(a):  // 내림차순
    a = cmp_swap(a, 0, 2)
    a = cmp_swap(a, 1, 3)
    a = cmp_swap(a, 0, 1)
    a = cmp_swap(a, 1, 2)
    return a

증명은  [1,2,3,4], [1,3,2,4], ....  [4, 3, 2, 1] ( 24개: 4*3*2*1 :  4! ) 같은 데이터 조합에 대해 모두 정렬이 된다면 잘 동작한다고 생각 할 수 있다. 좀 많이 무식하게 구현하면 프로스포츠 경기의 리그전 같이 하나의 element 가 모든 다른 element 와 비교하게 할 수도 있다. (아래에 position-counting-sort 이라고 정말 리그전 같은 알고리즘도 설명해 두었다.) 물론 sorting network 는 그 것보다는 좀 더 비교 연산이 적게 되어 있다. 

이 sorting network 를 예쁘게 만들면 아래처럼 된다. 

위 그림을 보면 4개의 비교 연산을 한 번에 하면 빠르게 동작 할 수 있겠다고 보여진다. 
1회   (a0,a1),  (a2,a3), (a4, a5), (a6, a7)
2회   (a0, a3), (a2, a3), (a4, a7), (a5, a6)
.....
6회  (a0,a1),  (a2,a3), (a4, a5), (a6, a7) 
로 계산할 수도 있다. 
그래픽 카드를 이용해서 병렬계산하거나 Simd(https://ko.wikipedia.org/wiki/SIMD) 명령어를 이용하면 고정적인 횟수로 빠르게 정렬을 할 수 있다. 

Timsort( https://en.wikipedia.org/wiki/Timsort )
* 기본 아이디어는 Insertion sort와 Merge sort를 결합한 알고리즘이다.  https://d2.naver.com/helloworld/0315536  여기에 너무 설명이 잘 되어 있다. 추가 설명이 필요 없을 것 같다.

pdqsort(https://github.com/orlp/pdqsort)
* golang 에서 이 sort 을 기본으로 사용 하겠다고 한다.  (https://github.com/golang/go/commit/72e77a7f41bbf45d466119444307fd3ae996e257 ) 
딱히 이 sort 에 대해 쉽게 설명이 잘 되어있는 곳을 찾지 못했다. 의도는 quick sort 와 quick sort 의 worst case 에 heap sort 를 하겠다는 발상이다. 


position-counting-sort(https://github.com/stgatilov/position-counting-sort) 
* 딱히 유명한 sort 방법은 아닌데, 요즘 개인적으로 관심을 가지는 방법이다. 위에서 설명한 Bitonic sorter 보다 더 무식하게 비교 연산을 한다.   a0, a1, a2, a3 이렇게 데이터가 있을 때,  a0 정렬 순위(정렬을 했을 때의 순위)계산하기 위해 싹 다 비교 연산을 한다. (마치 스포츠경기의 리그전 같이 한 팀이 나머지 팀을 모두 경기를 갖는 것이다. 거기에 따라서 순서가 정해진다. )

def position_counting_sort(a):  # 내림차순(높은 값에서 낮은 값으로)으로 정렬한다.
    ranks = [0] * len(a)
    for i in range(len(a)):
        ranking = 0
        for j in range(len(a)):
            # 무식하게 a[i] 보다 작은 것들의 개수를 센다. (구현이 귀찮아 자신도 포함해서 계산했다.)
            # 그 수 만큼 순위가 결정된다. 
            # a[i] 보다 큰 값이 하나도 없으면 ranking 은 0이다.
            if a[i] < a[j]:     
                ranking +=1
        # rank에 저장되는 값은 정렬된 순서가 아니라 각 위치별 순위이다.
        ranks[i] = ranking

    # rank 배열에 저장된 값에 따라 다시 한 번 여기서 순서대로 배열한다.
    out = [0] * len(a)
    for i in range(len(a)):
        out[ranks[i]] = a[i]
   return out

각 원소 하나에 대해 다른 원소들을 모두 비교한다. 일반적인 경우라면 당연히 매우 느릴 것이다. 그런데 simd 같이 한 번에 여러 비교 연산을 할 수 있다면 정렬 속도가 고정적인 괜찮은 알고리즘이 될 수 있다. 실용적인 관점에서 본다면 simd 같은 병렬처리를 사용하고 작은 그룹을 지어서 position-counting-sort 한 다음에 이 그룹끼리 merge 연산을 사용하는 hybrid 방식을 추구 해 볼 수 있다. 
물론 Bitonic sorter 보다 비교 연산이 많아서 불리 한 것 같지만 배열간의 이동이 거의 없다. 그리고 결과값으로 정렬된 데이터 뿐만 아니라 정렬순서표(위 소스에서 ranks 변수) 를 얻을 수 있다. 보통 key, value 가 하나의 set 로 정렬되는 경우도 많인데, 이 때 value 값에 정렬하면서 key 값을 얻을 수 있는 장점도 있다. 



 

   개인적으로 간단한 spreadsheet 기능을 가진 lightsheet_rust(https://github.com/yiunsr/lightsheet_rust ) 라는 프로그램을 개발중이다. 이 프로그램을 csv 파일을 열어서 sqlite3 로 저장해서 기본적인 sort, 수정, 추가, 간단한 수식을 제공하려고 한다.  이 프로그램 시작이 csv 파일을 열어서 sqlite3 DB에 저장하는 것이라보니 오래전 부터 https://yiunsr.tistory.com/821 이런 식으로 sqlite3 를 빠르게 insert 하는 것을 찾았다. 그러다 보니 이제는 rust 를 사용한다. 해당 option 외에 더 빠르게 하는 sqlite3 옵션은 없는 것 같다. 이제는 csv 파일을 빠르게 파싱하는 방법을 찾고 있다. 그러다 보니 simd (https://ko.wikipedia.org/wiki/SIMD) 를 이용하는 방법을 찾고 있다. 

 우선 빠른 csv 파서를 만들기 위해 필요한 것은 빠르게  field_separate(일반적으로 콤마, tsv 일 경우 탭), row_separate( \n, \r, \r\n) 그리고 특수한 문자 double quotes (큰따옴표 ") 를 빠르게 찾는게 필요하다. (큰따옴표가 의미 있는 이유는 a11, "b,11", c11 이라는 csv 파일 내용이 있으면 2번째 항목은 b,11  이다. 큰따옴표가 column 중간에 들어가는 콤마를 seperator 가 아닌 문자로 인식하게 해준다.)   일반적으로 파싱을 구현할 때 Tokenizer  로 token 을 분리할 때 문자 하나하나를 검사해서 구현한다. 그러나 csv 의 경우 특별한 의미를 가지는 위의 3가지 문자에 대해서만 빠르게 찾아서 Tokenizer 를 구현하게 된다면 빠르게 구현할 수 있다. 글로 설명하다보니 이해하기 어려워서 코드로 설명하면 

1) 문자하나하나 csv 를 파싱하는 방법

test_csv = "a1,b1,c1\na2,b2,c2"
for idx in range(len(test_csv)):
    ch = test_csv[idx]
    if ch == '"':
        quote 일 때 해야할 일
    elif ch == ",":
        comma 일 때 해야할 일
    elif ch == "\n":
        line feed 일 때 해야할 일
    else:
        pass
        # 아무 일도 안하고 위에서 많은 비교만 한다.

2) 특정 문자에 대해서 빠르게 이동 가능하다면

test_csv = "a1,b1,c1\na2,b2,c2"
for idx in find_3chars(test_csv, ",", '"',"\n"):
    ch = test_csv[idx]
    if ch == '"':
        quote 일 때 해야할 일
    elif ch == ",":
        comma 일 때 해야할 일
    elif ch == "\n":
        line feed 일 때 해야할 일
    else:
        pass
        # 여기가 불러질 일이 없다.

find_3chars 는 test_csv 에서 csv 를 파싱하는데 의미가 있는 문자 3가지( ",", '"',"\n") 의 index (여기서는  3,6, 9 ...)  를 찾아주는 함수이다. 

일반적으로 find_3chars 이런 함수는 일일이 문자를 찾는 것이라 2)가 1)과 같으면 같았지 더 빠를수가 없다. 그런데
intel simd 명령어중 avx2 에서만 동작하는 명령어를 사용하면 빠르게 find_3chars 를 구현 할 수 있다. simd 는 single instruction multiple data 의 약자다. 

 이 simd 명령어가 _mm256_cmpeq_epi8  (https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#text=_mm256_cmpeq_epi8&ig_expand=909 ) 이다. 이 명령어는  두 개의 256bit(1byte 32개) 피 연산자를 인자로 받아서 각 byte별로 결과를 256bit(32byte) 로 저장한다. 

a0  a1  a2 .... a31   (1byte 32개, 비교 대상으로 위에서의 a1,b1,c1\na2,b2,c2 )
b0  b1  b2 ... b31   (1byte 32개,     ,,,,, 가 32 개  )
c0  c1  c2  ... c31   (1byte 32개)

aN 과 bN 을 비교해서 같으면 cN은 0, 다르면 cN은 1 로 설정된다. 이 32개의 비교 연산이 하나의 instruction 으로 동작한다. 이 비교 연산을 하기 위해서는 b0 b1 ... b31 가 같은 값으로 32번 들어가야 한다. 
 a0 a1 a2 ... a31 내  "," 문자를 찾는다고 가정하면 b0 b1 b2 .... b31 모두 동일한 "," 가 들어가야 한다. 이 값을 넣는 것도 하나의 simd 명령어로 처리 가능하다.  _mm256_set1_epi8  라는 명렁어는 피연산자 one byte 를 받아서 32byte 에 동일하게 피연산자를 복사해서 넣는 처리를 해준다. 그리고  _mm256_cmpeq_epi8 를 통해 a0 ... a31 에 "," 문자가 있으면 c0 ... c31 중 해당하는 byte가  -1(0xFF)가 되고 아니면 0이 된다. 이 결과도 32byte가 아니라 결과를 bit 로 저장할 수 있다. _mm256_movemask_epi8 를 사용하면 32byte 의 결과를 다시 4byte(32bit) int 로 변경할 수 있다. ( _mm256_movemask_epi8 명령어는 msb(most significant bit) 가 있으면 1로 bit를 set 하고 없으면 0으로 bit 를 clear 한다. ).  이 결과를 대문자 C (32bit integer 이다. ) 라고 하자. 이렇게 하면  a0 .... a31 에서  콤마(,) 기호가 있으면 이 bit에 대응되는 C 에 bit가 1로 된다. 예를 들어 a1,b1,c12  => 2진수로 그 결과가     100100100  로 표시된다. 
 이렇게 \n 에 대해서도 찾고 '"' (큰따옴표)에 대해서도 찾아서 각각의 결과를  bitwise or  32byte 의 문자열에 대해 csv 파싱에 필요한 요소들의 위치를 한 번에 찾을 수 있다.  
 이렇게 찾은 방법이 simd 를 이용하기는 하지만 준비작업도 있고, 그 결과를 bit 로 된것을 분리해서 다시 확인하는 과정도 필요한다. 그래서 구현이 잘못되면 글자 하나하나를 비교하는 것보다 늦을 수도 있는데, 실제로 테스트 해보면 simd 를 이용한 방법이 확실히 더 빠르다. simd 를 이용한 방법이 뭔지 모르겠지만 메모리 cache 를 이용하면서 더 최적화되는 느낌이다. 

 내 경우 rust 를 이용해서 문자열에서 글자 3가지를 찾을 수 있는 라이브러리를 만들었다. https://github.com/yiunsr/bufchr  이다. 이 라이브러리는 https://github.com/BurntSushi/memchr 라는 라이브러리를
개선한 것이다. 기존의 memchr 이라는 라이브러리는 이미 계산된 32개의 결과를 모두 이용하지 않는다. 그냥 특정지점에서 부터 32byte 내에 매칭되는 문자가 있는지 확인하는 방식이다. 보통 찾는 곳은 haystack  이고 찾는 대상을 needle 이라고 표현한다. haystack 은 커다른 건초더미를 의미한다. 큰 건초더미에서 바늘을 찾는 것을 생각하면 된다. 기존 memchr 은 needle 이 적게 있을 때 빠른 속도를 낸다. 내가 만든 bufchr 은 바늘이 매우 많을 때 빠른 속도를 낸다. 특히 32byte 내에 needle 이 많으면 속도가 빠르게 된다. csv 의 경우 needle 이 32byte 내에 여러개 있을 가능성이 크기 때문에 내 라이브러리가 더 적합하다. 

 그리고 이 bufchr 를 이용해서 csv 파서를 만들고 있다. https://github.com/yiunsr/ss-csv 라는 rust 에서 csv 를 파서하는 라이브러리를 만들고 있다. 여기에는 cow string (clone on write 내용이 수정될 때에만 clone 됨)도 이용하고 있다. 따라서 기존https://github.com/BurntSushi/rust-csv 보다 더 빠른 속도를 보인다. 

테스트 해 보니 기존 것 보다 3분의 1 시간만에 csv 파싱을 끝내고 있다. 물론 csv 마다 다를 수 있다. 그리고 제약사항이 좀더 많다. 
https://gist.github.com/yiunsr/c0b0768d9e3938461214ec073f053b44   에서 테스트 코드를 확인 할 수 있다.