MonteCarloSimulation1

ConcurrentPerformBatch.swift

import Foundation

class ConcurrentPerformBatchExample {
    func isSatisfingF1(x: Double, y: Double) -> Bool {
        return x * x + (y - 5) * (y - 5) > 5 * 5
    }
    
    func isSatisfingF2(x: Double, y: Double) -> Bool {
        return (x * x) + (y * y) < 10 * 10
    }
    
    func isSatisfingF3(x: Double, y: Double) -> Bool {
        return (x - 5) * (x - 5) + (y - 10) * (y - 10) < 5 * 5
    }
    
    func isInsideS(x: Double, y: Double) -> Bool {
        return isSatisfingF1(x: x, y: y) && isSatisfingF2(x: x, y: y) && isSatisfingF3(x: x, y: y)
    }
    
    static let iterations = 10_000_0000
    var count = 0
    let dispatchQueueConcurrent = DispatchQueue.init(label: "worker1", attributes: .concurrent)
    let dispatchQueueSerial = DispatchQueue.init(label: "worker2")
    var sem = DispatchSemaphore(value: 1)
    
    func start() {
        let startTime = DispatchTime.now()
        DispatchQueue.concurrentPerform(iterations: ConcurrentPerformBatchExample.iterations) { [weak self] iter in
            guard let self = self else { return }
            let x = Double.random(in: 0.0...10.0)
            let y = Double.random(in: 0.0...10.0)
            
            if self.isInsideS(x: x, y: y) {
                _ = self.sem.wait(timeout: .distantFuture)
                self.count += 1
                self.sem.signal()
            }
            
            if iter % (ConcurrentPerformBatchExample.iterations/100) == 0 {
                print("\(iter / (ConcurrentPerformBatchExample.iterations/100)) / 100")
            }
        }
        let endTime = DispatchTime.now()
        let nanoSeconds = endTime.uptimeNanoseconds - startTime.uptimeNanoseconds
        let intervalTime = Double(nanoSeconds) / 1_000_000
        
        print("count = \(count)")
        print("S = \(Double(count) * 100 / Double(ConcurrentPerformBatchExample.iterations))")
        print("time = \(intervalTime) milsec")
        exit(EXIT_SUCCESS)
    }
}
/*
 iterations = 10_000_000
 count = 977229
 S = 9.77229
 time = 10743.776137
 
 iterations: 100_000_000
 count = 9776076
 S = 9.776076
 time = 148231.200488
 */

GPUComputationMTL.swift

import Foundation
import Metal
import MetalKit

class GPUComputationMTL {
    
    static let iterations = 100_000_000
    
    var startTime: DispatchTime?
    var endTime: DispatchTime?
    
    func start() {
        startTime = DispatchTime.now()
        guard let device = MTLCreateSystemDefaultDevice(),
              let defaultLibrary = device.makeDefaultLibrary(),
              let commandQueue = device.makeCommandQueue(),
              let isInsideSFunc = defaultLibrary.makeFunction(name: "isInsideS") else {
                print("Failed to create MTLdevice")
                return
        }
        var pipelineState: MTLComputePipelineState
        do {
            pipelineState = try device.makeComputePipelineState(function: isInsideSFunc)
        } catch {
            return
        }
        
        guard let commandBuffer = commandQueue.makeCommandBuffer(),
              let computeEncoder = commandBuffer.makeComputeCommandEncoder() else { return }
        
        var counter: Int32 = 0
        
        // setup MTLBuffer
        let bufferOut = device.makeBuffer(length: MemoryLayout<Bool>.size * GPUComputationMTL.iterations, options: .storageModeShared)
        let bufferCounter = device.makeBuffer(bytes: &counter, length: MemoryLayout<Int32>.size * 1, options: .storageModeShared)
        
        computeEncoder.setComputePipelineState(pipelineState)
        computeEncoder.setBuffer(bufferOut, offset: 0, index: 0)
        computeEncoder.setBuffer(bufferCounter, offset: 0, index: 1)
        
        let gridSize = MTLSizeMake(GPUComputationMTL.iterations / 1000, 1, 1)
        
        let threadGroupSize = MTLSizeMake(1000, 1, 1)
        computeEncoder.dispatchThreadgroups(gridSize, threadsPerThreadgroup: threadGroupSize)
        
        computeEncoder.endEncoding()
        commandBuffer.addCompletedHandler { _ in
            self.endTime = DispatchTime.now()
        }
        
        commandBuffer.commit()
        commandBuffer.waitUntilCompleted()

        
        guard let contents = bufferCounter?.contents() else { return }
        let ptr = contents.bindMemory(to: Int32.self, capacity: 1)
        let arr = UnsafeMutableBufferPointer(start: ptr, count: 1)


        var count = arr[0]
        
        print("count = \(count)")
        print("S = \(Double(count) * 100 / Double(GPUComputationMTL.iterations))")
        
        
        guard let startTime = startTime,
              let endTime = endTime else { return }
        
        let nanoSeconds = endTime.uptimeNanoseconds - startTime.uptimeNanoseconds
        let intervalTime = Float32(nanoSeconds) / 1_000_000
        print("time = \(intervalTime) milsec")
    }
}

/*
 GPUComputationMTL
 iterations = 100_000_000
 count = 9774663
 S = 9.774663
 time = 247.94894 milsec
 */

isInsideS.metal

//
//  isInsideS.metal
//  MonteCarloSimulation1
//
//  Created by KyungYoung Heo on 2021/09/11.
//

#include <metal_stdlib>
using namespace metal;

bool isSatisfingF1(float x, float y) {
    if (x * x + (y - 5) * (y - 5) > 5 * 5) {
        return true;
    } else {
        return false;
    }
}

bool isSatisfingF2(float x, float y) {
    if ((x * x) + (y * y) < 10 * 10) {
        return true;
    } else {
        return false;
    }
}

bool isSatisfingF3(float x, float y) {
    if((x - 5) * (x - 5) + (y - 10) * (y - 10) < 5 * 5) {
        return true;
    } else {
        return false;
    }
}

//float rand(int i){
//    return fract(sin(dot(i, float2(12.9898, 78.233))) * 43758.5453);
//}


uint rand_lcg(uint rng_state) {
    return 1664525 * rng_state + 1013904223;
}

uint rand_xorshift(uint rng_state) {
    // Xorshift algorithm from George Marsaglia's paper
    rng_state ^= (rng_state << 13);
    rng_state ^= (rng_state >> 17);
    rng_state ^= (rng_state << 5);
    return rng_state;
}

kernel void isInsideS(
            device bool* out,
            device atomic_uint* counter,
            uint index [[thread_position_in_grid]]
            ) {
    
    uint rng_state = rand_lcg(index);
    float x = float(rand_xorshift(rng_state)) * (1.0 / 4294967296.0) * 10;
    rng_state = rand_lcg(rng_state);
    float y = float(rand_xorshift(rng_state)) * (1.0 / 4294967296.0) * 10;
    
    
    out[index] = isSatisfingF1(x, y) && isSatisfingF2(x, y) && isSatisfingF3(x, y);
    if(out[index]) {
        atomic_fetch_add_explicit(counter, 1, memory_order_relaxed);
    }
}

main.swift

//
//  main.swift
//  MonteCarloSimulation1
//
//  Created by KyungYoung Heo on 2021/09/10.
//

import Foundation

//let semaphoreExample = SemaphoreExample()
//semaphoreExample.start()


//let concurrentPerformBatchExample = ConcurrentPerformBatchExample()
//concurrentPerformBatchExample.start()
//RunLoop.main.run()


let gPUComputationMTL = GPUComputationMTL()
gPUComputationMTL.start()

/*
 DispatchQueue.global().async
 count = 9768064
 S = 9.768064
 time = 967947.016009 milsec

 DispatchQueue.concurrentPerform
 count = 9776076
 S = 9.776076
 time = 148231.200488 milsec
 
 GPUComputationMTL
 count = 9774663
 S = 9.774663
 time = 247.94894 milsec
 */

SemaphoreExample.swift

import Foundation

class SemaphoreExample {
    func isSatisfingF1(x: Double, y: Double) -> Bool {
        return x * x + (y - 5) * (y - 5) > 5 * 5
    }

    func isSatisfingF2(x: Double, y: Double) -> Bool {
        return (x * x) + (y * y) < 10 * 10
    }

    func isSatisfingF3(x: Double, y: Double) -> Bool {
        return (x - 5) * (x - 5) + (y - 10) * (y - 10) < 5 * 5
    }

    func isInsideS(x: Double, y: Double) -> Bool {
        return isSatisfingF1(x: x, y: y) && isSatisfingF2(x: x, y: y) && isSatisfingF3(x: x, y: y)
    }

    static let iterations = 100_000_000
    var count = 0
    let dispatchQueueConcurrent = DispatchQueue.init(label: "worker1", attributes: .concurrent)
    let dispatchQueueSerial = DispatchQueue.init(label: "worker2")
    var sem = DispatchSemaphore(value: 1)
    let group = DispatchGroup()
    
    func start() {
        let startTime = DispatchTime.now()
        for i in 0..<SemaphoreExample.iterations {
            dispatchQueueConcurrent.async { [weak self] in
                guard let self = self else { return }
                let x = Double.random(in: 0.0...10.0)
                let y = Double.random(in: 0.0...10.0)

                if self.isInsideS(x: x, y: y) {
                    _ = self.sem.wait(timeout: .distantFuture)
                    self.count += 1
                    self.sem.signal()
        
                }
                
                if i % (SemaphoreExample.iterations / 100) == 0 {
                    print("\(i / (SemaphoreExample.iterations / 100)) / 100")
                }
            }
        }
        let endTime = DispatchTime.now()
        let nanoSeconds = endTime.uptimeNanoseconds - startTime.uptimeNanoseconds
        let intervalTime = Double(nanoSeconds) / 1_000_000
        
        print("count = \(count)")
        print("S = \(Double(count) * 100 / Double(SemaphoreExample.iterations))")
        print("time = \(intervalTime) milsec")
        exit(EXIT_SUCCESS)
    }
}
/*
 DispatchQueue.global().async
 iterations = 10_000_000
 count = 978294
 S = 9.78294
 time = 43135.978637
 
 iterations = 100_000_000
 count = 9768064
 S = 9.768064
 time = 967947.016009
 */

프로그래밍으로 이 문제를 풀어보면 어떄? 라는 얘기를 듣고 시작했다.

원래 문제와 풀이: https://www.youtube.com/watch?v=cPNdvdYn05c

이 문제를 푸는 두가지 방법

• 기하 해석
• 수치 해석

기하학적 방법: 부채꼴의 넓이에서 삼각형의 넓이를 빼는 방식으로 색칠된 넓이를 구할 수 있다.
수치해석 방법: 세 곡선의 함수를 구하고, 적분해서 색칠된 넓이를 구할 수 있다.

MonteCarlo Simulation으로 이 문제를 프로그래밍해서 풀 수 있다는 얘기를 듣고 주말에 코딩을 시작했다.
언어는 Swift, MacOS-Terminal App으로 프로젝트를 만들었다.

먼저 MonteCarlo Simulation하는 방법을 소개한다.

MonteCarlo 방법을 예를들어 소개하자면 사각형 안에 랜덤하게 돌을 던졌을 때 그 돌이 원 안에 들어갈 횟수를 세어서 원의 넓이의 근사값을 구하는 방법이다.
랜덤한 400개의 좌표를 생성하면 약 314개는 원 내부에 존재할 것이다.

이 시뮬레이션을 컴퓨터로 무수히 많이 반복할 수 있다.
랜덤한 (x, y) 좌표를 생성하고 원의 방정식을 이용해 좌표가 원 내부에 있는지 판별한다. 만약 원 내부에 있다면 count를 1 증가시킨다.
모든 반복이 끝나면 count를 총 반복횟수 N으로 나눈다.

친구가 파이썬으로 구한 원의 넓이 근사값

원래의 문제를 좌표계위에 올려놓으면 다음과 같이 나타낼 수 있다.

세 곡선으로 둘러싸인 영역 S위의 한 점 (x, y)는 원의 방정식을 이용해 나타낼 수 있다.
f1의 바깥쪽 f2, f3의 안쪽에 포함되는 점은 수식으로 표현된다.

이 문제는 이제 한변의 길이가 10인 정사각형 안에서 랜덤한 좌표를 생성할 때 S의 영역 안에 좌표가 몇 번 포함되는지 계산하는 문제로 치환할 수 있다.

[0,10] 범위의 랜덤한 (x, y)좌표를 생성하고, isSatisfyingF1(x, y) && isSatisfyingF2(x, y) && isSatisfyingF3(x, y)가 true이면 그 좌표는 S의 영역 안에 포함되므로 count를 증가시킨다. 이 과정을 1억 번 반복하고 S = count * 100 / 1억 을 계산하면 영역 S의 넓이의 근사값을 구할 수 있다.

main.swift에서는 이 값을 세가지 방법을 이용해 계산한다.

SemaphoreExample.swift는 DispatchQueue와 DispatchSemaphore를 이용해 1억번 반복하면서 시뮬레이션 하는 코드이다.
Swift의 GCD에서는 block된 쓰레드 수가 증가할수록 성능이 저하되는데 (https://developer.apple.com/documentation/dispatch/dispatch_queue의 Avoiding Excessive Thread Creation를 참고) 쓰레드가 30개까지 생성되고 갑자기 모든 쓰레드가 blocked 되면서 한동안 아무 작업도 하지 않는다. 한참 뒤에는 몇개의 쓰레드가 종료되고 남은 쓰레드에서 작업을 실행하다가 다시 쓰레드가 일정 개수 생성되면 멈추는 것을 반복한다.
DispatchQueue.global()를 써서 쓰레드 성능을 높일 수 있다는 코멘트도 있어서 시도해봤지만 큰 차이는 없었다. 글로벌 큐나 직접 생성한 큐나 같은 현상이 발생했다. 그래서 1억번 반복했을 때는 약 967초가 걸렸다.

이 현상을 개선하기 위해 여러 방법을 시도해보고 개선한 코드가 ConcurrentPerformBatch.swift이다.
DispatchQueue.concurrentPerform(iterations:)을 사용할 때는 쓰레드가 8~10개만 처음에 생성되고 추가로 생기거나 사라지지 않았다.
내 생각엔 쓰레드풀에 계속 작업을 넣어주는 오버헤드가 감소해서 전체적인 성능이 올라간 것 같다.
1억번 반복했을 때 시간도 약 148초로 큰 개선이 있었다.

속도가 빨라지는게 재미있어서 조금 더 개선을 해보기로 했다. 이 코드는 부동소수점 연산이 많고 꽤 많은 반복작업을 하기 때문에 GPU연산을 이용해보기로 했다. GPU연산은 이런 작업을 하기에 최적화 되어있다. 코드는 GPUComputationMTL.swift과 isInsideS.metal를 참고하기 바란다.

Apple은 Metal이라는 그래픽 프레임워크를 사용한다. GPUComputationMTL.swift에서Compute Shader 셋업을 하고 isInsideS.metal에서는 Shader language로 이 문제를 구현했다.
이 과정에서 어려웠던 점은 1억개의 결과 중 true에 해당하는 횟수만큼 counter를 증가시키는 법이었다. 병렬로 실행되는 쓰레드에서 한 변수를 atomic하게 acc하는 방법을 찾느라 헤매었다. atomic_fetch_add_explicit() 를 이용해 함수 이름처럼 atomic하게 fetch하고 add해서 정확한 횟수를 구할 수 있었다.
그 다음으로는 Metal에서 random한 좌표를 생성하는 부분이었다. OpenGL에서 사용하는 코드를 참고해서 gid를 시드넘버로 이용해 pseudo random한 좌표를 생성해냈다.

처음엔 이 문제를 흥미로 접근했다가 Swift에서 구현하는 과정에서 여러가지 문제를 경험하고 생각보다 오래 걸려서 끝냈다. 더 많은 개선점이 남아있는 코드지만 일단 여기에서 마친다.

그래서 이 문제의 원래 답은 다음 그림과 같고, arctan값을 구해서 계산해보면 색칠된 영역의 넓이는 9.773570749999998가 나온다.

1억번 반복했을때 영역 안에 포함된 횟수 count
프로그래밍해서 계산한 넓이 S와 걸린시간은 다음과 같다.

 DispatchQueue.global().async
 count = 9768064
 S = 9.768064
 time = 967947.016009 milsec

 DispatchQueue.concurrentPerform
 count = 9776076
 S = 9.776076
 time = 148231.200488 milsec
 
 GPUComputationMTL
 count = 9774663
 S = 9.774663
 time = 247.94894 milsec

랜덤하게 좌표를 생성하기 때문에 실행할 때마다 결과가 조금씩 달라지긴 하지만 대체로 비슷한 값이 나온다.
걸린 시간이 967초 -> 148초 -> 0.24초로 개선된 건 주목할만한 것 같다.