0

I am trying to make a histogram of Poisson random generated variables using Python and C. I wanted to use Python for plotting and C for generating. This resulted in the following to codes

Python:

import ctypes
import numpy as np
import matplotlib.pyplot as plt
import time

lam = 5.0
n = 1000000

def generate_poisson(lam, n):
    array = np.zeros(n, dtype= np.int)
    f = ctypes.CDLL('./generate_poisson.so').gen_poisson
    f(ctypes.c_double(lam), ctypes.c_int(n), ctypes.c_void_p(array.ctypes.data))
    return array

start_time = time.time()
array = generate_poisson(lam,n)
print(time.time() - start_time)

plt.hist(array, bins = [0,1,2,3,4,5,6,7,8,9,10,11,12], density = True)
plt.savefig('fig.png')
print(array)

C:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

double get_random() { return ((double)rand() / (double)RAND_MAX); }

int poisson_random(double lam){
    int X;
    double prod, U, explam;

    explam = exp(-lam);
    X = 0;
    prod = 1.0;
    while (1){
        U = get_random();
        prod *= U;
        if (prod > explam){
            X+=1;
        }
        else {
            return X;
        }
    }
}

void gen_poisson(double lam, int n, void * arrayv)
{
    int * array = (int *) arrayv;
    int index = 0;
    srand(time(NULL));

    for (int i =0; i<n; i++, index++){
        //printf("before %d\n", array[i]);
        array[index++] = poisson_random(lam);
        //printf("after %d\n", array[i]);
    }
}

The problem understanding why this works, or at least it seems to work correctly, occurs inside the for loop in gen_poisson(). Somehow using array[index++] instead of array[index] results in a correct histogram. But I dont really understand why this has to be the case. The code also seems to work when the for loop is changed to

for (int i =0; i<2*n; i++){
        //printf("before %d\n", array[i]);
        array[i++] = poisson_random(lam);
        //printf("after %d\n", array[i]);
}

Can someone explain why in this case the loop has to be incremented twice? I just started programming in C while I have some experience in Python. So assume the culprit is my lack of understanding about C. Thank you in advance

10
  • Can you clarify in how far it does not work without post-incrementing index? The title and tags mention Cython, which generally makes writing correct code much easier, but the question does not seem to include it – do you actually use Cython? In how far is it relevant for your problem? Commented Oct 25, 2020 at 13:39
  • You've got an array of C long, not int. Commented Oct 25, 2020 at 13:40
  • Also, NumPy already has a Poisson distribution sampling function. Commented Oct 25, 2020 at 13:42
  • I actually meant ctypes instead of Cython. I edited the post, my bad Commented Oct 25, 2020 at 13:44
  • 1
    Your Python code declared an np.array of c_int...so changing to long* in the C code makes no sense to me. Your code worked as is for me, with int array and int* with the correct indexing. Commented Oct 26, 2020 at 0:26

1 Answer 1

1

Changing the gen_poisson into:

void gen_poisson(double lam, int n, void * arrayv)
{
    long * array = (long *) arrayv;
    srand(time(NULL));
    for (int i =0; i<n; i++){
        array[i] = poisson_random(lam);
    }
}

solves the problem. The problem was as pointed out with declaring the array as int * instead of long *.

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.