/*
* Alex Laubscher
* Gillespie Algorithm
* Uses a GPU generator for the numbers
*/

#include <curand.h>
#include <stdio.h>
#include <time.h>

int main() {

    // Initializing variables for gillespie algorithm
    int counter;
    int death;
    int total;
    double tau;
    double sample;
    int check;

    // Initialize variables for the GPU generator
    int count = 2500000;
    curandGenerator_t gen;
    float *devURN;
    float *hostURN;

    // Allocate n floats on host
    hostURN = (float *)calloc(count, sizeof(float));

    // Allocate n floats on device
    cudaMalloc((void **) &devURN, count*sizeof(float));

    // Create the generator
    curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT);

    // Set the seed
    curandSetPseudoRandomGeneratorSeed(gen, 1234ULL);

    // Initial population
    int pop = 0;

    // Initializing time
    double time = 0;
    double maxTime = 100000;

    // Can be outside cuz it never changes
    int birth = 1000;

    // Starting the timer
    clock_t time_elapsed = clock();

    // Run the while loop over 100,000 simulation seconds
    while (time < maxTime) {
        // Setting the propensity of the rxn
        death = pop;

        // Sum over the propensities
        total = birth + death;

        // Need to cast the double
        check = counter % (count / 2);

        if (check == 0) {
            // Generate the floats
            curandGenerateUniform(gen, devURN, count);

            // Copy the numbers back to the device
            cudaMemcpy(hostURN, devURN, count*sizeof(float),
                cudaMemcpyDeviceToHost);
        }

        // Calculate time step
        tau = (1.0 / total) * log(hostURN[check * 2]);

        // Second random choice
        sample = total * (hostURN[check * 2 + 1]);

        // Update populations based on second urn
        if (sample < birth) {
            pop = pop + 1;
        } else {
            pop = pop - 1;
        }

        // Update the time step
        time = time - tau;

        // Increment the counter
        counter++;
    }

    // End the time and convert to sec
    time_elapsed = clock() - time_elapsed;
    double timer = ((double) time_elapsed) / CLOCKS_PER_SEC;

    //Calculate the reactions per sec
    double rate = counter / timer;
    printf("Population: %f\n", pop);
    printf("Counter: %d\n", counter);
    printf("Timer: %f\n", timer);
    printf("Rate: %f\n", rate);

    curandDestroyGenerator(gen);
    cudaFree(devURN);
    free(hostURN);

    return 0;
}

	.file	"tmpxft_002f7cd6_00000000-6_gillespieGPU.cudafe1.cpp"
	.text
#APP

#NO_APP
	.type	_ZL26__cudaUnregisterBinaryUtilv, @function
_ZL26__cudaUnregisterBinaryUtilv:
.LFB2030:
	.cfi_startproc
	endbr64
	movq	_ZL20__cudaFatCubinHandle(%rip), %rdi
	jmp	__cudaUnregisterFatBinary@PLT
	.cfi_endproc
.LFE2030:
	.size	_ZL26__cudaUnregisterBinaryUtilv, .-_ZL26__cudaUnregisterBinaryUtilv
	.section	.rodata.str1.1,"aMS",@progbits,1
.LC5:
	.string	"Population: %f\n"
.LC6:
	.string	"Counter: %d\n"
.LC7:
	.string	"Timer: %f\n"
.LC8:
	.string	"Rate: %f\n"
	.section	.text.startup,"ax",@progbits
	.globl	main
	.type	main, @function
main:
.LFB2027:
	.cfi_startproc
	endbr64
	pushq	%r15
	.cfi_def_cfa_offset 16
	.cfi_offset 15, -16
	movl	$4, %esi
	movl	$2500000, %edi
	pushq	%r14
	.cfi_def_cfa_offset 24
	.cfi_offset 14, -24
	pushq	%r13
	.cfi_def_cfa_offset 32
	.cfi_offset 13, -32
	pushq	%r12
	.cfi_def_cfa_offset 40
	.cfi_offset 12, -40
	xorl	%r12d, %r12d
	pushq	%rbp
	.cfi_def_cfa_offset 48
	.cfi_offset 6, -48
	pushq	%rbx
	.cfi_def_cfa_offset 56
	.cfi_offset 3, -56
	subq	$56, %rsp
	.cfi_def_cfa_offset 112
	movq	%fs:40, %rax
	movq	%rax, 40(%rsp)
	xorl	%eax, %eax
	call	calloc@PLT
	leaq	32(%rsp), %rdi
	movl	$10000000, %esi
	movq	%rax, %rbp
	call	cudaMalloc@PLT
	leaq	24(%rsp), %rdi
	movl	$100, %esi
	call	curandCreateGenerator@PLT
	movq	24(%rsp), %rdi
	movl	$1234, %esi
	call	curandSetPseudoRandomGeneratorSeed@PLT
	call	clock@PLT
	movq	$0x000000000, 8(%rsp)
	movq	%rax, %r14
.L7:
	movl	%r13d, %eax
	movl	$1250000, %ecx
	leal	1000(%r12), %r15d
	cltd
	idivl	%ecx
	movl	%edx, %ebx
	testl	%edx, %edx
	jne	.L3
	movq	32(%rsp), %rsi
	movq	24(%rsp), %rdi
	movl	$2500000, %edx
	call	curandGenerateUniform@PLT
	movq	32(%rsp), %rsi
	movl	$2, %ecx
	movq	%rbp, %rdi
	movl	$10000000, %edx
	call	cudaMemcpy@PLT
.L3:
	addl	%ebx, %ebx
	movslq	%ebx, %rbx
	movss	0(%rbp,%rbx,4), %xmm0
	call	logf@PLT
	cvtsi2sdl	%r15d, %xmm2
	movsd	.LC1(%rip), %xmm1
	cvtss2sd	%xmm0, %xmm0
	divsd	%xmm2, %xmm1
	movss	.LC2(%rip), %xmm2
	mulsd	%xmm0, %xmm1
	cvtsi2ssl	%r15d, %xmm0
	mulss	4(%rbp,%rbx,4), %xmm0
	comiss	%xmm0, %xmm2
	jbe	.L11
	incl	%r12d
	jmp	.L6
.L11:
	decl	%r12d
.L6:
	movsd	8(%rsp), %xmm3
	movsd	.LC3(%rip), %xmm0
	incl	%r13d
	subsd	%xmm1, %xmm3
	comisd	%xmm3, %xmm0
	movsd	%xmm3, 8(%rsp)
	ja	.L7
	call	clock@PLT
	cvtsi2sdl	%r13d, %xmm1
	movl	%r12d, %edx
	leaq	.LC5(%rip), %rsi
	subq	%r14, %rax
	movl	$2, %edi
	cvtsi2sdq	%rax, %xmm0
	divsd	.LC4(%rip), %xmm0
	xorl	%eax, %eax
	movsd	%xmm0, 8(%rsp)
	divsd	%xmm0, %xmm1
	movq	%xmm1, %rbx
	call	__printf_chk@PLT
	movl	%r13d, %edx
	movl	$2, %edi
	xorl	%eax, %eax
	leaq	.LC6(%rip), %rsi
	call	__printf_chk@PLT
	movsd	8(%rsp), %xmm0
	leaq	.LC7(%rip), %rsi
	movb	$1, %al
	movl	$2, %edi
	call	__printf_chk@PLT
	movq	%rbx, %xmm0
	leaq	.LC8(%rip), %rsi
	movb	$1, %al
	movl	$2, %edi
	call	__printf_chk@PLT
	movq	24(%rsp), %rdi
	call	curandDestroyGenerator@PLT
	movq	32(%rsp), %rdi
	call	cudaFree@PLT
	movq	%rbp, %rdi
	call	free@PLT
	movq	40(%rsp), %rax
	subq	%fs:40, %rax
	je	.L8
	call	__stack_chk_fail@PLT
.L8:
	addq	$56, %rsp
	.cfi_def_cfa_offset 56
	xorl	%eax, %eax
	popq	%rbx
	.cfi_def_cfa_offset 48
	popq	%rbp
	.cfi_def_cfa_offset 40
	popq	%r12
	.cfi_def_cfa_offset 32
	popq	%r13
	.cfi_def_cfa_offset 24
	popq	%r14
	.cfi_def_cfa_offset 16
	popq	%r15
	.cfi_def_cfa_offset 8
	ret
	.cfi_endproc
.LFE2027:
	.size	main, .-main
	.type	_ZL24__sti____cudaRegisterAllv, @function
_ZL24__sti____cudaRegisterAllv:
.LFB2053:
	.cfi_startproc
	endbr64
	pushq	%rax
	.cfi_def_cfa_offset 16
	leaq	_ZL15__fatDeviceText(%rip), %rdi
	call	__cudaRegisterFatBinary@PLT
	movq	%rax, %rdi
	movq	%rax, _ZL20__cudaFatCubinHandle(%rip)
	call	__cudaRegisterFatBinaryEnd@PLT
	leaq	_ZL26__cudaUnregisterBinaryUtilv(%rip), %rdi
	popq	%rdx
	.cfi_def_cfa_offset 8
	jmp	atexit@PLT
	.cfi_endproc
.LFE2053:
	.size	_ZL24__sti____cudaRegisterAllv, .-_ZL24__sti____cudaRegisterAllv
	.section	.init_array,"aw"
	.align 8
	.quad	_ZL24__sti____cudaRegisterAllv
	.section	.nvFatBinSegment,"aw"
	.align 8
	.type	_ZL15__fatDeviceText, @object
	.size	_ZL15__fatDeviceText, 24
_ZL15__fatDeviceText:
	.long	1180844977
	.long	1
	.quad	fatbinData
	.quad	0
	.local	_ZL20__cudaFatCubinHandle
	.comm	_ZL20__cudaFatCubinHandle,8,8
	.section	.rodata.cst8,"aM",@progbits,8
	.align 8
.LC1:
	.long	0
	.long	1072693248
	.section	.rodata.cst4,"aM",@progbits,4
	.align 4
.LC2:
	.long	1148846080
	.section	.rodata.cst8
	.align 8
.LC3:
	.long	0
	.long	1090021888
	.align 8
.LC4:
	.long	0
	.long	1093567616
	.ident	"GCC: (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0"
	.section	.note.GNU-stack,"",@progbits
	.section	.note.gnu.property,"a"
	.align 8
	.long	1f - 0f
	.long	4f - 1f
	.long	5
0:
	.string	"GNU"
1:
	.align 8
	.long	0xc0000002
	.long	3f - 2f
2:
	.long	0x3
3:
	.align 8
4:

/*
* Alex Laubscher
* Gillespie Algorithm
* Uses a GPU generator for the numbers
*/


#include <hip/hip_runtime.h>
#include <hiprand/hiprand.h>
#include <stdio.h>
#include <time.h>

int main() {

    // Initializing variables for gillespie algorithm
    int counter;
    int death;
    int total;
    double tau;
    double sample;
    int check;

    // Initialize variables for the GPU generator
    int count = 2500000;
    hiprandGenerator_t gen;
    float *devURN;
    float *hostURN;

    // Allocate n floats on host
    hostURN = (float *)calloc(count, sizeof(float));

    // Allocate n floats on device
    hipMalloc((void **) &devURN, count*sizeof(float));

    // Create the generator
    hiprandCreateGenerator(&gen, HIPRAND_RNG_PSEUDO_DEFAULT);

    // Set the seed
    hiprandSetPseudoRandomGeneratorSeed(gen, 1234ULL);

    // Initial population
    int pop = 0;

    // Initializing time
    double time = 0;
    double maxTime = 100000;

    // Can be outside cuz it never changes
    int birth = 1000;

    // Starting the timer
    clock_t time_elapsed = clock();

    // Run the while loop over 100,000 simulation seconds
    while (time < maxTime) {
        // Setting the propensity of the rxn
        death = pop;

        // Sum over the propensities
        total = birth + death;

        // Need to cast the double
        check = counter % (count / 2);

        if (check == 0) {
            // Generate the floats
            hiprandGenerateUniform(gen, devURN, count);

            // Copy the numbers back to the device
            hipMemcpy(hostURN, devURN, count*sizeof(float),
                hipMemcpyDeviceToHost);
        }

        // Calculate time step
        tau = (1.0 / total) * log(hostURN[check * 2]);

        // Second random choice
        sample = total * (hostURN[check * 2 + 1]);

        // Update populations based on second urn
        if (sample < birth) {
            pop = pop + 1;
        } else {
            pop = pop - 1;
        }

        // Update the time step
        time = time - tau;

        // Increment the counter
        counter++;
    }

    // End the time and convert to sec
    time_elapsed = clock() - time_elapsed;
    double timer = ((double) time_elapsed) / CLOCKS_PER_SEC;

    //Calculate the reactions per sec
    double rate = counter / timer;
    printf("Population: %f\n", pop);
    printf("Counter: %d\n", counter);
    printf("Timer: %f\n", timer);
    printf("Rate: %f\n", rate);

    hiprandDestroyGenerator(gen);
    hipFree(devURN);
    free(hostURN);

    return 0;
}

	.text
	.file	"gillespieGPU.hip"
	.section	.rodata.cst8,"aM",@progbits,8
	.p2align	3, 0x0                          # -- Begin function main
.LCPI0_0:
	.quad	0xbff0000000000000              # double -1
.LCPI0_2:
	.quad	0x40f86a0000000000              # double 1.0E+5
.LCPI0_3:
	.quad	0x412e848000000000              # double 1.0E+6
	.section	.rodata.cst4,"aM",@progbits,4
	.p2align	2, 0x0
.LCPI0_1:
	.long	0x447a0000                      # float 1000
	.text
	.globl	main
	.type	main,@function
main:                                   # @main
	.cfi_startproc
# %bb.0:
	pushq	%rbp
	.cfi_def_cfa_offset 16
	pushq	%r15
	.cfi_def_cfa_offset 24
	pushq	%r14
	.cfi_def_cfa_offset 32
	pushq	%r13
	.cfi_def_cfa_offset 40
	pushq	%r12
	.cfi_def_cfa_offset 48
	pushq	%rbx
	.cfi_def_cfa_offset 56
	subq	$40, %rsp
	.cfi_def_cfa_offset 96
	.cfi_offset %rbx, -56
	.cfi_offset %r12, -48
	.cfi_offset %r13, -40
	.cfi_offset %r14, -32
	.cfi_offset %r15, -24
	.cfi_offset %rbp, -16
	movl	$2500000, %edi                  # imm = 0x2625A0
	movl	$4, %esi
	callq	calloc
	movq	%rax, %rbx
	leaq	24(%rsp), %rdi
	movl	$10000000, %esi                 # imm = 0x989680
	callq	hipMalloc
	leaq	32(%rsp), %r14
	movq	%r14, %rdi
	movl	$400, %esi                      # imm = 0x190
	callq	hiprandCreateGenerator
	movq	(%r14), %rdi
	movl	$1234, %esi                     # imm = 0x4D2
	callq	hiprandSetPseudoRandomGeneratorSeed
	xorl	%r14d, %r14d
	callq	clock
	xorpd	%xmm0, %xmm0
	movq	%rax, %r15
                                        # implicit-def: $ebp
.LBB0_1:                                # =>This Inner Loop Header: Depth=1
	movsd	%xmm0, 16(%rsp)                 # 8-byte Spill
	movslq	%ebp, %r12
	imulq	$1801439851, %r12, %rax         # imm = 0x6B5FCA6B
	movq	%rax, %rcx
	shrq	$63, %rcx
	sarq	$51, %rax
	addl	%ecx, %eax
	imull	$1250000, %eax, %eax            # imm = 0x1312D0
	subl	%eax, %r12d
	jne	.LBB0_3
# %bb.2:                                #   in Loop: Header=BB0_1 Depth=1
	movq	32(%rsp), %rdi
	movq	24(%rsp), %rsi
	movl	$2500000, %edx                  # imm = 0x2625A0
	callq	hiprandGenerateUniform
	movq	24(%rsp), %rsi
	movl	$10000000, %edx                 # imm = 0x989680
	movq	%rbx, %rdi
	movl	$2, %ecx
	callq	hipMemcpy
.LBB0_3:                                #   in Loop: Header=BB0_1 Depth=1
	leal	1000(%r14), %r13d
	xorps	%xmm0, %xmm0
	cvtsi2sd	%r13d, %xmm0
	movsd	.LCPI0_0(%rip), %xmm1           # xmm1 = mem[0],zero
	divsd	%xmm0, %xmm1
	movsd	%xmm1, 8(%rsp)                  # 8-byte Spill
	leal	(%r12,%r12), %eax
	cltq
	xorps	%xmm0, %xmm0
	cvtss2sd	(%rbx,%rax,4), %xmm0
	callq	log
	mulsd	8(%rsp), %xmm0                  # 8-byte Folded Reload
	movsd	16(%rsp), %xmm2                 # 8-byte Reload
                                        # xmm2 = mem[0],zero
	addsd	%xmm0, %xmm2
	xorps	%xmm0, %xmm0
	cvtsi2ss	%r13d, %xmm0
	leal	1(,%r12,2), %eax
	cltq
	mulss	(%rbx,%rax,4), %xmm0
	xorl	%eax, %eax
	movss	.LCPI0_1(%rip), %xmm1           # xmm1 = mem[0],zero,zero,zero
	ucomiss	%xmm0, %xmm1
	seta	%al
	leal	(%r14,%rax,2), %r14d
	decl	%r14d
	incl	%ebp
	movsd	.LCPI0_2(%rip), %xmm0           # xmm0 = mem[0],zero
	ucomisd	%xmm2, %xmm0
	movapd	%xmm2, %xmm0
	ja	.LBB0_1
# %bb.4:
	callq	clock
	subq	%r15, %rax
	xorps	%xmm1, %xmm1
	cvtsi2sd	%rax, %xmm1
	divsd	.LCPI0_3(%rip), %xmm1
	movsd	%xmm1, 8(%rsp)                  # 8-byte Spill
	xorps	%xmm0, %xmm0
	cvtsi2sd	%ebp, %xmm0
	divsd	%xmm1, %xmm0
	movsd	%xmm0, 16(%rsp)                 # 8-byte Spill
	movl	$.L.str, %edi
	movl	%r14d, %esi
	xorl	%eax, %eax
	callq	printf
	movl	$.L.str.1, %edi
	movl	%ebp, %esi
	xorl	%eax, %eax
	callq	printf
	movl	$.L.str.2, %edi
	movsd	8(%rsp), %xmm0                  # 8-byte Reload
                                        # xmm0 = mem[0],zero
	movb	$1, %al
	callq	printf
	movl	$.L.str.3, %edi
	movsd	16(%rsp), %xmm0                 # 8-byte Reload
                                        # xmm0 = mem[0],zero
	movb	$1, %al
	callq	printf
	movq	32(%rsp), %rdi
	callq	hiprandDestroyGenerator
	movq	24(%rsp), %rdi
	callq	hipFree
	movq	%rbx, %rdi
	callq	free
	xorl	%eax, %eax
	addq	$40, %rsp
	.cfi_def_cfa_offset 56
	popq	%rbx
	.cfi_def_cfa_offset 48
	popq	%r12
	.cfi_def_cfa_offset 40
	popq	%r13
	.cfi_def_cfa_offset 32
	popq	%r14
	.cfi_def_cfa_offset 24
	popq	%r15
	.cfi_def_cfa_offset 16
	popq	%rbp
	.cfi_def_cfa_offset 8
	retq
.Lfunc_end0:
	.size	main, .Lfunc_end0-main
	.cfi_endproc
                                        # -- End function
	.type	.L.str,@object                  # @.str
	.section	.rodata.str1.1,"aMS",@progbits,1
.L.str:
	.asciz	"Population: %f\n"
	.size	.L.str, 16

	.type	.L.str.1,@object                # @.str.1
.L.str.1:
	.asciz	"Counter: %d\n"
	.size	.L.str.1, 13

	.type	.L.str.2,@object                # @.str.2
.L.str.2:
	.asciz	"Timer: %f\n"
	.size	.L.str.2, 11

	.type	.L.str.3,@object                # @.str.3
.L.str.3:
	.asciz	"Rate: %f\n"
	.size	.L.str.3, 10

	.type	__hip_cuid_,@object             # @__hip_cuid_
	.bss
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               # 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_

	.text
	.p2alignl 7, 3214868480
	.fill 96, 4, 3214868480
	.type	__hip_cuid_,@object             ; @__hip_cuid_
	.section	.bss,"aw",@nobits
	.globl	__hip_cuid_
__hip_cuid_:
	.byte	0                               ; 0x0
	.size	__hip_cuid_, 1

	.ident	"AMD clang version 18.0.0git (https://github.com/RadeonOpenCompute/llvm-project roc-6.3.2 25012 e5bf7e55c91490b07c49d8960fa7983d864936c4)"
	.section	".note.GNU-stack","",@progbits
	.addrsig
	.addrsig_sym __hip_cuid_
	.amdgcn_target "amdgcn-amd-amdhsa--gfx1100"
	.amdgpu_metadata
---
amdhsa.kernels:  []
amdhsa.target:   amdgcn-amd-amdhsa--gfx1100
amdhsa.version:
  - 1
  - 2
...

	.end_amdgpu_metadata
