Adding very large arrays

Question

I am writing a stochastic simulation for a Yule process:

If you start with a certain number of bins each containing a random number of balls, you add another ball to an existing bin with a probability that is proportional to the number of balls within that bin, as well as create another bin with a single ball each time step.

Because I am trying to basically weed out any stochastic noise, I am using very large number of iterations, on the order of 1e6. This code was originally developed in MATLAB, making much less complicated, but is far too slow for large simulations.

Is there a better way of doing this?

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

void init_KISS();
unsigned int JKISS();
double uni_dblflt();
double array_sum();

int main()
{
    int t,count,i,j,k,ii,jj,kk;
    int desval;
    desval=1000000;
    unsigned int r;
    double v,u,bsum;
    init_KISS();    
    for(i=0; i<100000;i++)
    {
        r = JKISS(); /* waste the first 100000 random numbers */
        uni_dblflt();
    }
    /*allocate memory for all params*/
    double *bins=malloc((desval+1)*sizeof(double));/*Urns, keeps track of actual "balls" in each bin*/
    if(bins==NULL)
    {
        printf("out of memory\n");
        exit(1);
    }

    double *p=malloc((desval+1)*sizeof(double));/* Urns/(sum of all "balls) i.e. the probability
                                                 that a new ball get placed in the ith urn. */
    if(p==NULL)
    {
        printf("out of memory\n");
        exit(1);
    }

    double *ptot=malloc((desval+1)*sizeof(double));/*Sum of all probabilities for the 1st to the ith
                                                    bin*/
    if(ptot==NULL)
    {
        printf("out of memory\n");
        exit(1);
    }

    for(i=0;i<(desval+1);i++)
    {
        bins[i]=0;
        ptot[i]=0;
        p[i]=0;
    }

    for(i=0;i<10;i++)/*seed 10 initial bins with 1 "ball"*/
    {
        v=1;
        bins[i]=v;
    }


    for(j=10;j<desval;j++)
    {
        bsum=array_sum(bins,desval);/*sum the balls in each bin*/
        for(k=0;k<j;k++) /* calculate the probability that with which a ball might be placed in 
                          each urn*/
        {
            p[k]=bins[k]/bsum;
        }

        for(ii=0;ii<desval;ii++)/*Evaluate the probability total cutoffs*/
        {
            ptot[ii]=array_sum(p,ii);
        }
        u=uni_dblflt();
        count=0;
        while(u<=ptot[count])/*choose bin to place new ball. The numerical value of ptot(count) 
                              is directly related to the probability space with which u can be 
                              picked*/
        {
            count++;
        }

        bins[count]++;
        bins[j]=1;
    }

    for(jj=0;jj<desval;jj++)
    {
        printf("%.8f\n",bins[jj]);
    }

    free(bins);
    free(p);
    free(ptot);
    bins=NULL;
    p=NULL;
    ptot=NULL;
    return 0;
}

double array_sum(double *arr, int maxval)
{
    /* this function sums a large array*/
    int i;
    double sum; 
    sum=0;

    for(i=0; i<maxval;i++)
    {
        sum +=arr[i];
    }

    return sum;
}

Answer 1

This is really an algorithms question. For example, consider this loop:

for(ii=0;ii<desval;ii++)/*Evaluate the probability total cutoffs*/
{
    ptot[ii]=array_sum(p,ii);
}

why are you calling array_sum() every time through this loop when the only thing that changes is the addition of one more element? Why can't this loop be:

sum= 0;
ptot[0]= 0;
for(ii=1;ii<desval;ii++)/*Evaluate the probability total cutoffs*/
{
    ptot[ii]=(sum+= p[ii - 1]);
}

Answer 2

Let's get some preliminary problems out of the way before we discuss efficiency — and there is much efficiency to be gained.

Random number generator

You appear to be using the advice and public-domain code from David Jones for your pseudo-random number generator. Basically,

• init_KISS() seeds the pseudorandom number generator
• JKISS() generates a random unsigned int
• uni_dblflt() generates a random double between 0 and 1 using a uniform distribution.

It would have been nice to note that information in the question.

Declutter your main() by moving this code out into a function with a self-documenting name.

void initialize_prng()
{
    init_KISS();
    for (int discard = 0; discard < 100000; discard++)
    {
        JKISS(); /* waste the first 100000 random numbers */
        uni_dblflt();
    }
}

Proliferation of variables

First, let the compiler alert you to problems. Compiling with warnings turned on, we see problems about how variables are being used.

$ gcc -Wall -c yule.c
yule.c:13:9: warning: unused variable 't' [-Wunused-variable]
    int t,count,i,j,k,ii,jj,kk;
        ^
yule.c:13:29: warning: unused variable 'kk' [-Wunused-variable]
    int t,count,i,j,k,ii,jj,kk;
                            ^
2 warnings generated.

For that matter, r is assigned but never used, so it could also be eliminated. v is just a constant 1, and can also be eliminated.

The overall problem here is that you have a proliferation of variables that makes the code hard to follow. Ever since C99, you can declare variables anywhere inside a function, not just at the top. (It's now 2014 — surely it's OK to use C99?) Best practice is to declare your variables as close to the point of use as possible, as it improves readability, decreases mental workload, and puts variables in a tighter scope to avoid silly mistakes.

Other matters of style

An int is only guaranteed to hold 16 bits, or numbers up to 32767. In practice, your code will work because an int will be 32 bits on any modern machine. However, if you want to be certain that you can store 10⁶, strictly speaking, you should use a long.

Think about readability when using loop-counting variables. For example, for(j=10;j<desval;j++) is notionally a continuation of for(i=0;i<10;i++), just with a different loop body. Therefore, I would use i as the counter for the second loop as well.

I don't see any need here for #include <math.h> and #include <time.h>.

Instead of zeroing the arrays manually, use calloc(), which zeroes out the memory it allocates.

You allocate one more element than you need for each array, as if making a zero-terminator element at the end. That's something that one would do with NUL-terminated strings and some arrays of pointers. For this problem, a zero-terminator is useless, since 0.0 is a legal value that can't act as a sentinel. Therefore, the code would be less puzzling if you dropped the pretense of having a zero-terminator.

Declutter your main() function further by extracting the printing loop into its own function.

Efficiency

As @llewelly noted, you are doing unnecessary work to sum the arrays repeatedly. How bad is it?

for(j=10;j<desval;j++)
{
    …

    for(ii=0;ii<desval;ii++)/*Evaluate the probability total cutoffs*/
    {
        ptot[ii]=array_sum(p,ii);
    }
    …
}

Three nested loops makes it O(n³). When n = 10⁶, that's significant!

@llewelly suggests removing the innermost array_sum(), to make it O(n²). We can actually do much, much better than that.

First, note that the p array is a useless intermediary between bins and ptot.

Next, consider combining bins and ptot into a single array. Suppose that instead of using bins to store the number of balls in each element, we had a cumulative_bins that stored the cumulative number of balls from bin 0 to bin i. Reconstructing what bins would have held is trivial: just subtract the predecessor from cumulative_bins[i]. Reconstructing what ptot would have held is also easy: divide cumulative_bins[i] by cumulative_bins[desval - 1], where cumulative_bins[desval - 1] is the total number of balls in the system. (Note that this cumulative_bins[] should be discrete, not floating-point.)

If you try to implement that idea, you'll find that reading the cumulative number of balls in bins 0 to i is an O(1) operation, adding one ball to a bin is an O(n) operation, and finding the index of a cumulative count is O(n). Fortunately, there is a nice data structure that can improve on that: a binary indexed tree, also known as a Fenwick Tree, which can do all three operations in O(log n) time.

---

Solution

This program performs the job in O(n log n) time, in contrast to your original O(n³). Excluding the time it takes to print the output, it completes a simulation of n = 10⁶ bins in about half a second.

long array_get_cumulative(const long fenwick_tree[], long idx)
{
    long sum = 0;
    for (; idx >= 0; idx = (idx & (idx + 1)) - 1)
    {
        sum += fenwick_tree[idx];
    }
    return sum;
}

void array_incr(long fenwick_tree[], long idx, size_t size, long incr)
{
    for (; idx < size; idx |= (idx + 1))
    {
        fenwick_tree[idx] += incr;
    }
}

long array_find_cumulative(long fenwick_tree[], long cumFre, size_t size)
{
    long idx = 0, tIdx, bitMask;
    for (bitMask = 1; size & ~((bitMask << 1) - 1); bitMask <<= 1) {}
    for (; bitMask && ((tIdx = idx + bitMask) < size); bitMask >>= 1)
    {
        if (cumFre >= fenwick_tree[tIdx - 1])
        { 
            // if current cumulative frequency is equal to cumFre, 
            // we are still looking for higher index (if exists)
            idx = tIdx;
            cumFre -= fenwick_tree[tIdx - 1];
        }
    }
    return idx - (cumFre == 0);
}

void array_print(long fenwick_tree[], size_t size)
{
    for (long i = 0, preceding_cumulative = 0; i < size; i++)
    {
        long cumulative = array_get_cumulative(fenwick_tree, i);
        printf("%ld\n", cumulative - preceding_cumulative);
        preceding_cumulative = cumulative;
    }
}

int main()
{
    long n = 1000000L;

    initialize_prng();

    /* Urns, keeps track of cumulative "balls" in each bin (including all preceding bins) */
    long *bins = calloc(n, sizeof(long));
    if (bins == NULL)
    {
        fprintf(stderr, "Out of memory\n");
        exit(1);
    }

    for (long i = 0; i < 10; i++) /* seed 10 initial bins with 1 "ball" each */
    {
        array_incr(bins, i, n, 1);
    }

    for (long i = 10; i < n; i++)
    {
        long ball_count = array_get_cumulative(bins, n - 1);
        double u = uni_dblflt();

        long j = array_find_cumulative(bins, u * ball_count, n);

        array_incr(bins, j, n, 1);
        array_incr(bins, i, n, 1);
    }

    array_print(bins, n);
    free(bins);
    return 0;
}

There may be off-by-one errors in array_find_cumulative() — I'll leave it to you to work them out.