An Intersection of Politics and Optimization

I had a thought a few weeks ago, that something today jogged my memory about, and that is that one of the problems that most people agree on with regards to politics is gerrymandering. I've seen a few ways to try to combat it, from getting an independent (if you can find one) panel to divvy up a state, to just creating square blocks within a state and calling that good.

I have my own idea, but first I need to set it up so that it makes sense. As I see it, there are a number of competing, and possibly mutually exclusive, objectives when creating districts, and that is before even thinking about the political persuasion of the constituents. These objectives include:

• Equal population in each district
• District lines fall on county boundaries
• Non-gerrymandered
• Human Input into final decision (nobody likes to be told what to do by a computer)

There may be others, but I think these are a good start. It is almost guaranteed that optimizing for one will necessarily not be optimized for the others. As it is, we are currently optimized only on the last point, which we really can't assign an optimization function to anyway.

The first three, however, I think we can assign optimization values to. The first one is pretty easy, for a state of population P with N districts, where the i^th district has a population p_i then minimize the value of Sum(abs(p_i - P/N)). You can use a squared value instead of abs if you want, that's just using the L2 norm instead of the L1 norm. Anyway, minimize that sum, and populations in each district will be equalized. For population data, use census data by precinct, assume that the population is uniform across the precinct (inaccurate, I know, but a good first guess).

The second one is a bit harder, but I think one way to express it mathematically is to minimize the the area which can be drawn between a district boundary and a county line. For a county that is split in such a way, count only the smaller area (otherwise the county could be cut any way you like). I'm sure others could come up with a way to optimize this as well that might make a bit more sense.

The third one could be mathematically imposed by minimizing the sum of the differences between the convex hull for a district and the district itself. If all of the districts are convex, then this sum will be 0. You will also have a non-gerrymandered state. For states such as Maryland, with its jagged southern border, it would probably be impossible to get this sum to 0, but that's still OK, we just want to minimize it for any given state.

Now, as I said before, these three cannot be complete optimized at the same time, and different people will put different weights on each of the objectives. Some people might believe that the population equality is most important, while others might believe that districts following county lines is paramount. This is almost where the human element comes in. See, even though people might disagree on what is most important, what can be done before assigning weights is to try to optimize over all of them. What will result is something called a Pareto frontier. Each point on the Pareto frontier represents a different weighting scheme. Some reasonable limits might also want to be enforced, maybe that every p_i has to fall within 0.75 * P/N to 1.25 * P/N.

What would then happen is that a multiple objective optimization program can find a bunch of points on the Pareto frontier, in this case, each point is actually a specific map of the state in question. A set of these can then be handed off to a state legislature, or some other deciding body, and then they can select the one they like the most, bringing the human element back into it. They can then select the most optimal map for them if they like.

So there it is, a possible way to reduce gerrymandering while still leaving the final decision in the hands of people, possibly even state legislatures. I even think that this is a non-partisan idea. Also, I'd be happy to hear other (non-partisan) objectives that such a program might want to optimize for. It would be interesting if multiple states did this to see what the relative weights each state chose for it's given district mapping.

Goals

So yeah, this thing that I never update.

Interesting thing about last year's posts, besides the fact that I didn't make very many of them. The most popular one was Building gcc 4.8. Mostly because it actually got linked in a Stack Overflow answer about the same thing. That was a point of pride for me :).

Looking forward to the future though, I have a few goals. Posting more on here is not one of them. I may post more, I may not. Last year I kind of stopped posting because I moved in July, and then I moved again in November. Between packing, unpacking, and unpacking again (movers did the packing for the second move), this kind of got left by the wayside. Anyway, on to some of my actual goals for 2014:

Fitness:

• Do 2000 pushups, this actually seems pretty low. But if I hit it, I can always up the count next year
• Spend 2000 minutes on exercise machines probably elipticals, though I think I'll count an actual bicycle toward this as well
• Run a 5K, no walking. My eventual goal is a half-marathon (and maybe a full marathon someday, but we'll start with the half). But I want achievable goals for 2014, and I'm not sure the half-marathon is achievable for me

Household:

• Build a 6 foot tall bookshelf. I built a solid pine 4' x 4' bookshelf a couple years ago. It is the strongest, and best looking bookshelf we own. It costs a lot more than the cheapo particle board ones, but the end result, and the experience, are worth it.
• Play 36 board games during the year with my son. It is a way for me to spend quality time with my child as well as building his love for a hobby I love.

Technology/Professional:

• Get a MySQL certification, I have to do more and more with databases, getting a certification is probably a good idea, at least from a resume perspective.
• Submit at least one conference or journal paper where I am the lead author. I've got one in mind, but I'm not sure it is publication worthy
• Complete at least 1 online course and get the certificate, I'm signed up for Convex Optimization through Stanford online, you can join me there if you'd like.

Well, there they are. They should be achievable. I'll start on them tomorrow, or this weekend at the latest :)

CUDA - Figuring out how many threads have launched Part 2

Part one is here.

Last time I had determined that I could see how many threads had been launched on my GTS 450.  However, each one of those threads was making an atomicAdd() call to a pointer in host memory. Two things that  are quite bad.

If I want the host to be able to read what the threads are doing, then I need to keep the memory location in host.  However, I don't want to have every single thread updating that memory location. Thankfully, I don't have to.  Blocks are loaded onto individual SMPs at once and retired all at once. Not only that, but I can get the block size quite easily from inside the kernel with the blockDim variable.  Implementing this change, my new kernel looks like this:

    __global__ void wasteTime(int nIter, int *threads){

        int blkThreads = 0;
        if(     threadIdx.x == 0 &&
                threadIdx.y == 0 &&
                threadIdx.z == 0){
            blkThreads = blockDim.x * blockDim.y * blockDim.z;
            atomicAdd(threads, blkThreads);
        }

        int j;
        float data = 1.0;
        for (j = 0; j < nIter; j++){
            data = cos(data) + data;
        }

        if(     threadIdx.x == 0 &&
                threadIdx.y == 0 &&
                threadIdx.z == 0){
            atomicSub(threads, blkThreads);
        }
    }

I also made a modification to the main routine after the kernel call to be able to make more updates as to how many threads have been launched on the device:

        wasteTime<<<grid, block, 0, stream>>>(nIter, threads);

        struct timespec delay;
        delay.tv_sec = 0;
        delay.tv_nsec = 200000000L;

        printf("%d threads running\n", *threads);
        nanosleep(&delay, NULL);

        while (*threads > 0){
            printf("%d threads running\n", *threads);
            nanosleep(&delay, NULL);
        }

        cudaFree(threads);

This also requires a #include <time.h> up at the top of the source. One last change is that I swapped the grid to be 4x4x4, to launch more blocks.  The full code is here. Running the code generated the following results:


0 threads running
2048 threads running
2048 threads running
2048 threads running
1024 threads running

As you can see, the first time I called the printf statement, no threads hadn't been spawned.  The total number of threads launched was actually (4 x 4 x 4) x (4 x 4 x 4) = 4096 , but  I apparently got a maximum number of threads (for this problem and GPU) that could be resident on the GPU as being 2048.  Another good thing, however, was that I can see the number of threads decreasing at the right right before it finishes.

There is a small problem in the kernel, especially with small block sizes though.  Any time there is a divergence in threads within a warp (set of 32 threads), all of the threads in that warp perform all of the commands for all of the divergent paths with the non-applicable threads having their results masked out.  What this means for my code is that threads 1-31 will all wait and twiddle their thumbs while thread 0 goes off and does its DMA atomicAdd and atomicSub. With the size of my block equaling 64, this means that a full half of my threads are affected by the time delay associated with 1/64 of my threads.

I can attack this in a few ways. The first is to do the atomic commands only on a single thread from the entire grid. Another is to set the value of *threads prior to the launch and remove the atomicAdd, but leave in the atomicSub. I think I want to leave the atomicSub at the end of the kernel anyway (and calculate blkThreads within that if statement) because I like seeing individual blocks retired. Changing the kernel to have only block 0,0,0 thread 0,0,0 atomicAdd the number of threads for the grid and then have thread 0,0,0 for each block atomicSub it back down (code here) has the following results:

0 threads running

4096 threads running

3072 threads running

2048 threads running

1024 threads running

The count starts off at 4096 and goes down as the blocks retire. Having the host increment the number of threads prior to launch, though, may be even better.  This would allow me to see how many threads have been scheduled for launch, which is what I really want to know for load balancing purposes.

To avoid race conditions I'll have to still use atomicAdd, just from the host, especially if I branch out into threading, where we may have multiple threads trying to do this all at once. Simple right? Well, you can try it, but I'll spoil it, atomicAdd is a "device" function, it can only be called from the device. Obviously the idea of atomics aren't new to CUDA though, so we can look else where.

I could try to use the C++11 standard library atomics, but there are a few  but that requires me to create the variables as a std::atomic<int> object. This is problematic because I also need to allocate the memory using cudaMallocHost to get my DMA access.  Hmm, well, I can also use intrinsics. GCC has an intrinsic __sync_fetch_and_add.  This works in C (which means I can drop it into my runningThreads.cu file), but it is, as I said, intrinsic to GCC, making it less portable. This could bite me later if I want to use the Intel or PGI compilers.

I implemented the increment with __sync_fetch_and_add and changed the loop so that it would launch a new grid (with a total of 3 grids being launched) when there was less than 30% of the threadsPerGrid still on the device. This could just as easily be set based on some percentage of cudaDeviceProp::maxThreadsPerMultiProcessor * cudaDeviceProp::multiProcessorCount. In any case, the code can be found here, and the output was:

16384 threads running

16384 threads running

12288 threads running

12288 threads running

8192 threads running

8192 threads running

20480 threads running

20480 threads running

16384 threads running

16384 threads running

12288 threads running

12288 threads running

8192 threads running

8192 threads running

20480 threads running

20480 threads running

16384 threads running

16384 threads running

12288 threads running

12288 threads running

8192 threads running

8192 threads running

4096 threads running

4096 threads running

Looking pretty good. Next time I'll throw some OpenMP threads into the mix and maybe an idea I just had about a way I could use the C++11 atomics to capture the number of currently launched threads, but I haven't even tried to implement that yet.

CUDA - Figuring out how many threads have launched Part 1

I've talked about multiple threads sharing a single GPU once before. That got me thinking though, what if I don't want all of them to be using the GPU simultaneously.  What if the GPU actually becomes the bottleneck in my computational code? In all likelyhood, anything I implement on a GPU will probably also have a CPU implementation somewhere, because I would like my code to run on any system. If it can run on an accelerated system, all the better.

Now that I've laid out my motivation, it's time to figure out what I can do about it.  It is pretty simple to figure out how much of the devices memory has been utilized with the function cudaMemGetInfo(). But what about the number of blocks or threads? Now, it may be that my Google-Fu is weak, but I couldn't find a built-in way to let me determine the current number of threads running, so I decided to see if I could build one myself.

I started off using the same code I used last time, though I modified my  wasteTime routine so that it would cut back on the registers used a bit.  The first iteration of the new wasteTime function looks like this:

     __global__ void wasteTime(int nIter, int *threads){

        atomicAdd(threads, 1);

        int j;
        float data = 1.0;
        for (j = 0; j < nIter; j++){
            data = cos(data) + data;
        }
    }

In the main function, I got rid of the OpenMP (for now) to focus on getting this call right.

     int main(int argc, char *argv[]){
        int nIter;
        nIter = 100000;
        int *threads;
        cudaMallocHost(&threads, sizeof(int));

        *threads = 0;

        dim3 grid(2, 1, 1);
        dim3 block(4, 4, 4);

        cudaStream_t stream;

        cudaStreamCreate(&stream);

        wasteTime<<<grid, block, 0, stream>>>(nIter, threads);

        printf("%d threads started\n", *threads);
        cudaStreamDestroy(stream);
        printf("%d threads started\n", *threads);

        cudaFree(threads);
    }

A few of notes before I run this:

1. The threads pointer in main is the same threads pointer in the CUDA code, this is done via pinned memory on the host and direct memory access (DMA) from the GPU to the host.
2. DMA from the device to the host is incredibly expensive.
3. DMA is not available for all compute cards, though it is for all Fermi-based cards and above, and honestly, Fermi is old enough now that you should be able to assume that if people are going to accelerate their code, they'll have it.
4. Atomic operations are expensive because they get serialized.
5. I have changed the grid and block dimensions for launch so that they can be a bit more general. I get the feeling I'll be using this general project as a framework for a bunch of different tests. If I don't, no biggie.
6. The source code up to this point is here.

I've left in the streams from last time because I want to parallelize this later.

Building and running (still on my 4 core AMD with a GTS 450) I get an answer of:

0 threads started
71 threads started

Hmm, interesting. Doing some basic math real quick, we know that we have 2 blocks, and each block has 4*4*4 = 64 threads. The first time that the printf statement is called, we haven't actually launched any threads yet, or at least they haven't modified the global memory yet. We then destroy the stream, which I had thought would synchronize the device and the host, and when we call printf again, we only have 71 threads, more than 1 block, but not fully 2.

Now we'll go back and synchronize the device to the host. This will be a simple change where we just add a call to cudaDeviceSynchronize() right after the second printf call, and we'll put another printf right after that to see how things have changed. Also, I've modified the wasteTime function and given it an atomicAdd(threads, -1); at the end, so that instead of seeing how many threads have started, I can see how many threads are currently running.  I've also modified the printf statements to reflect this. The code for the new version is here. The results from 3 separate runs are:

ross@kronix:~/dev/cuda/runningThreads$ Debug/runningThreads 
0 threads running
113 threads running
0 threads running
ross@kronix:~/dev/cuda/runningThreads$ Debug/runningThreads 
0 threads running
77 threads running
0 threads running
ross@kronix:~/dev/cuda/runningThreads$ Debug/runningThreads 
0 threads running
68 threads running
0 threads running

As you can see there is a fair amount of entropy in the number of threads "caught" running by the middle printf, which is between the thread destroy and the device synchronize.

I definitely have some more to add to this, such as not having every thread run the atomic functions, and spawning into multiple OpenMP threads as well. I'll be following this post up with some future posts detailing my thought process and approach to solving this problem I may or may not have invented out of thin air.

In fact, you may be wondering why I even care about this at all. The reason is I may be working on a project soon which will literally have thousands of concurrent tasks each one of which may or may not be able to saturate a GPU for a short time. If I am even running on an accelerated system, it will still have a high CPU to GPU ratio on each node. While spawning some of these tasks onto the GPU to free up CPU time will be beneficial, if I end up being stuck waiting for the GPU to crunch through all of the tasks I've given it, because I haven't paid attention to load balancing, well, that would be the opposite of accelerated functionality, now wouldn't it. That last sentence is long and has a lot of commas, but I'm going to leave it like that anyway.

On a more pragmatic tangent, like "why wouldn't you just develop for CPUs and then accelerate with a Xeon Phi?", The fact that I can do these kinds of tests and experiments at home is one of the nice things about how NVIDIA has set up their system. I can tinker with my low cost system, but port that code and expertise over to the production systems at work with minimum effort. While I'd love to get my hands on an Intel Xeon Phi and do some work with it (a minimum $2k proposition at the low end, including the fact that my current desktop couldn't run it), I would pretty much have to purchase the Intel Compilers to go with it to be able to build toward it. Contrast that with my $65 video card and free compilers, and the choice isn't that hard at the moment. That and I find it fun.

Still, let me reiterate that: I would love to get my hands on an Intel Xeon Phi.

Finally, on a non-GPU related note, I want to thank Tod for giving me some space on his server. He set up my account a while ago for D&D stuff, but I don't think he'll mind if I put some code up for others to try out. If you want me to take it down (or actually link to you, other than tod.net) let me know.

Thoughts on Spring Cleaning

I feel like my wife and I were pretty productive today. First I finally got around to mowing the lawn for the first time this season. There was one point where the grass was so thick and I was moving too fast (at a walk) that I actually stalled the mower. Going over that patch a bit slower let me take it out though.

After that, and lunch, we started in on what we refer to as the office. This is where our computers are and it was in pretty bad shape. Truth be told it still isn't actually clean, but I can actually see the surfaces of both of the desks that make up my corner. I also raised my two monitors (having experienced 2 monitors, I far prefer the increased real estate to a single monitor) up onto a mostly unused shelf that was sitting behind them.

This made an incredible difference for my entire desk. Not only did it make the space they were sitting in usable, it made the mostly unused space behind them (and under the shelf) usable as well. Now I'll just have to get used to the raised position of the monitors. While I am slouching pretty bad right now, for the most part I have found that I have had better posture with them standing higher.

I'm sure this post has been pretty boring so far, so here are some amusing tidbits from our adventures in cleaning.

1. We found a small LEGO creator kit, I built the race car and it is now sitting on the shelf with my monitors, I may or may not have been playing with the race car all evening to the annoyance of my wife.
2. We filled a medium sized trash bin twice with recyclables.
3. We found a birthday card that was supposed to go to somebody a few weeks ago (technically this was found in the kitchen, and they did get the present at the appropriate time).
4. We found a hospital bill that was over a year old.
5. Pants which had a loose button, and therefore were on the sewing table for over year were fixed.
6. ???
7. Profit

Now to get ready for the D&D adventure I'm running tonight.

Sharing a (nVidia) GPU between OpenMP Threads

In case you couldn't tell from the title, this will be somewhat code heavy.

I was thinking about one of my problems for work and thought that at some point it may be worth while for me to know whether or not I could share a GPU between different threads using OpenMP. I thought it was possible, since they would all share the same context (still not sure what that means, but I   know that MPI processes do not share the same context, and that they can't share a GPGPU for pre-Kepler hardware).

A couple of disclaimers, I only write in CUDA when I'm doing GPGPU work, which I know ties me to nVidia hardware, but that's the way it is.  Also, this code has no error checking, because it was a quick and dirty test this evening. I'm pretty sure I spent more time on writing this blog post than I did on  writing the code.

At work I have access to most compilers as well as a GTX 680 (in my workstation, it's sweet), some K10s, a K20, a system where I can get up to 4 nodes with each node having one M2050 (Fermi), and a system with multiple M2070s per node (Fermi). At home though, well, I'll let some of the device properties speak for themselves.

 
==== Device 0 ===
                  name: GeForce GTS 450
    concurrent kernels: T
                   SMP: 4
   Max Threads per SMP: 1536
  Max Theads per Block: 1024
            Global Mem: 1.00 G
               Compute: 2.1
   Asyncronous Engines: 1

So yeah, not the most powerful GPU out there, but it has computer capacity 2.1, which was why I bought it. It only cost ~$65 at the time. But I think my next one might be in the GTX Titan family. But that is a different post.

As an aside, if anybody wants the code to generate the above queries the code is here.

The first thing that I decided I had to do was create a kernel which could waste some serious time so that I wouldn't have the possibility that one kernel from a thread would be done before the overhead from the next even started. I didn't want to deal with allocating and freeing memory, though. I came up with a very simple and useless kernel:

 
    __global__ void wasteTime(int nIter){
        int j;
        float data = 1.0;
        for (j = 0; j < nIter; j++){
            data = cos(data) + sin(data);
        }
    }

Next, I set up a simple program to call the kernel from a bunch of different threads.

 
    int main(int argc, char *argv[]){
        int nIter;
        nIter = 100000;

    #pragma omp parallel
        {
            wasteTime<<<1,128>>>(nIter);
        }
    }

I created this project using the nVidia Nsight Eclipse Edition (I'm running on the Linux box), which has nvvp (a handy-dandy GPGPU profiler) built in, which will allow me to easily see if I was able to overlap my kernel calls. Here is the a profile of that first attempt

Figure 1.

So the good news is that my kernel takes more than half a second, which makes it a supreme time waster for this situation.  Also, at first glance, it appears as though I have 3 kernels that are overlapping. First glance is wrong though, and my initial inkling that it is wrong was that I have 4 cores on my CPU, so I should see 4 kernels if everything was working correctly. Looking a bit closer, I'm pretty sure I actually have only 1 kernel running there. To check I put in a printf statement inside the parallel section. If I do have more than 1 thread running, I should see more than one output line.

    int main(int argc, char *argv[]){
        int nIter;
        nIter = 100000;

    #pragma omp parallel
        {
            wasteTime<<<1,128>>>(nIter);
            printf("done.\n");
        }
    }

Figure b.

Hmm, only one 'done.' That's a problem. And then I remember, I have to pass in some flags to make sure that it actually compiles the #pragma omp parallel statement. I didn't see any warnings about not being able to handle the #pragma statements, which I usually get. That may have something to do with nvcc, or it may have something to do with me only checking to see if successfully built, and not looking through the console output.  It took me a bit to figure out where they go when using Nsight Eclipse edition (I often role my own Makefiles), but I'll give you a couple of screen shots so that you can easily find them if you do the same thing.  You can get the dialog below by selecting 'Project' ->'Properties'.

Figure iii.

Note, I only opened the dialog once, I did not have 2 dialogs one on top of the other, just thought I would make sure that was clear. I then had to go into my Debug folder and call 'make clean' (though there was probably another way to do this), rebuild and reprofile, yielding the following results:

Figure we're getting closer.

Well, we got the 'done' to show up 4 times, and I can also see the 4 individual threads in the profiler as well.  The bad news is that there is absolutely no overlap of the time wasting kernels.  There should be plenty of room on the GPU, after all, there are 4 SMP, and each one can handle 1.5k threads. I've also established that it can handle concurrent kernels.

And that's when I remember it, streams.  At least, I think that is the answer. Basically all of the threads are calling into the GPU's default stream (which appears to be stream 2 according to the profiler). Not only that, the default stream doesn't do things quite as asynchronously as other streams. To be able to get concurrent kernels I have to have each of the threads call into a different stream. That should be easy enough.

    int main(int argc, char *argv[]){
        int nIter;
        nIter = 100000;

        cudaStream_t stream;

    #pragma omp parallel
        {
            cudaStreamCreate(&stream);
            wasteTime<<<1,128, 0, stream>>>(nIter);
            cudaStreamDestroy(stream);
            printf("done.\n");
        }
    }

Figure Whoo popcorn superposition (kernel overlap)!

I've collapsed the CPU threads as they aren't all that important for what I want to test. However, as we can see stream 7 and stream 10 are both executing at the same time. The red line in there is apparently overhead directly associated with the profiler, and I may look into that some more.  I also ran the (from the debug directory) /urs/bin/time -p ./ompTest which gave me an elapsed time of 1.07 seconds. That is obviously more than 1x the time to run wasteTime(), but less than 2x the time to run wasteTime() so the profiler might be messing up some other stuff in there too.

Finally, final version of my (somewhat small) source code file is here. If you've made it this far, hopefully this post has been of some use to you.

Thoughts on exercise

I don't get to exercise as much as I would like to. Or rather, I don't make as much time for exercise as I aught to. I enjoy exercise in general. In the past I have enjoyed soccer, ultimate, and even a bit of basketball (though I cannot claim to ever have had much skill in that last one). For solo exercise, I actually somewhat enjoy using an elliptical, and I think that I would like to run at least a half-marathon at some point. Maybe I should make myself an official 'bucket list', but that is not what this post is about.

The primary problems for me are motivation and time management. In terms of motivation, I almost always enjoy the exercise while I am doing it, and I typically feel pretty good afterward as well. For some reason, though, it can be hard to drag myself out of bed or away from the computer to actually start the activity.

The issue of time management has become more pronounced lately. I used to hit the gym on Monday, Wednesday, and Friday mornings. However, about a month ago we got a dog. What this means, though, is that I typically take him on a walk every morning and evening, which takes about 20 minutes. So while I do get some exercise every day, it is not very intense. I have considered running a bit on these walks, and I may start doing that.

Another option is to start going to the gym in the evenings, but that takes time away from my family. So this post comes to an end without any real resolution, but rather just some thoughts which have been bouncing around my head lately.