Calculating PI in C# using Monte-Carlo simulation

The following code sample shows numeric compuation of the number PI using Monte-Carlo simulation. First, a sequential approach is used. Then the Parallel.For construct in TPL is used. In the end, we use Tasks in TPL.


To compute PI we use the same approach. We consider the unit circle inscribed in a square around origo with corners at coordinates (-1,-1), (-1, 1), (1,1) and (1,-1). The number PI can be defined as generating random numbers and looking at the ratio of the numbers inside the circle M divided upon the total numbers generated N. We know that the following can then be expected:

(a) M / N = PI / 4

Why? Because the square has got an area equal to four, remember that the unit square got sides equal to the number 2 and its area is therefore 2 * 2 = 4. The unit circle got a radius of 1, hence its area is PI * 1^2 = PI. The ratio to be expected between the areas of the unit circle and unit rectangle therefore gives the formula above. We can further compute the approximated numeric value of PI equal to:

(b) PI = 4 * (M / N)

This expression (b) is directly from the previous expression (a)
Let's move on the code sample, review the code. I have included a screen shot at the end. The conclusion I got after testing showed after several runs shows that the sequential version runs in about 3.5 seconds on my eight core system with about half the time, about 1.8 seconds using Parallel.For - The last version using Tasks and Tasks.WaitAll give about 1.7 seconds and the quickest compuation, about twice as fast. The iterations I used in the demo was 80 million.
Here is the code written in C#:

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace MonteCarloPiApproximation
{
    class Program
    {
        static int numberOfCores = Environment.ProcessorCount; 
        static int iterations = 10000000 * numberOfCores;

        static void Main(string[] args)
        {
            Console.WriteLine("Monte Carlo numeric simulation of PI");
            Console.WriteLine("Iteration limit: " + iterations);
            Console.WriteLine("Number of processor cores on system: " + Environment.ProcessorCount);
            var sw = new Stopwatch();
            sw.Start();

            Console.WriteLine("\nMONTE CARLO SIMULATION");

            MonteCarloPiApproximationSerialSimulation();
            sw.Stop();

            Console.WriteLine("Serial simulation: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Restart();

            Console.WriteLine("\nMONTE CARLO SIMULATION");
            MonteCarloPiApproximationParallellForSimulation(); 

            sw.Stop();

            Console.WriteLine("Parallell simulation using Parallel.For: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Restart();

            Console.WriteLine("\nMONTE CARLO SIMULATION");

            MonteCarloPiApproximationParallelTasksSimulation(); 

            Console.WriteLine("Parallell simulation using parallell Tasks: (ms)" + sw.ElapsedMilliseconds);
            Console.WriteLine();

            sw.Stop();           

            Console.WriteLine("Press Enter Key");



            Console.ReadKey(); 
        }

        private static void MonteCarloPiApproximationParallelTasksSimulation()
        {
            double piApproximation = 0;
            int inCircle = 0;
            double x, y = 0;

            int[] localCounters = new int[numberOfCores];
            Task[] tasks = new Task[numberOfCores];

            for (int i = 0; i < numberOfCores; i++)
            {
                int procIndex = i; //closure capture 
                tasks[procIndex] = Task.Factory.StartNew(() =>
                {
                    int localCounterInside = 0;

                    Random rnd = new Random();

                    for (int j = 0; j < iterations / numberOfCores; j++)
                    {
                        x = rnd.NextDouble();
                        y = rnd.NextDouble();
                        if (Math.Sqrt(x * x + y * y) <= 1.0)
                            localCounterInside++;
                    } 
                    localCounters[procIndex] = localCounterInside;

                });               
            }

            Task.WaitAll(tasks);
            inCircle = localCounters.Sum(); 

            piApproximation = 4 * ((double)inCircle / (double)iterations);

            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));
           
        }      

        private static void MonteCarloPiApproximationParallellForSimulation()
        {
            double piApproximation = 0;
            int inCircle = 0;
            double x, y = 0;
                   
            Parallel.For(0, numberOfCores, new ParallelOptions{ MaxDegreeOfParallelism = numberOfCores }, i =>
            {
              
                int localCounterInside = 0;

                Random rnd = new Random(); 

                for (int j = 0; j < iterations / numberOfCores; j++)
                {
                    x = rnd.NextDouble();
                    y = rnd.NextDouble();
                    if (Math.Sqrt(x*x+y*y) <= 1.0)
                        localCounterInside++;                                                        
                }

                Interlocked.Add(ref inCircle, localCounterInside); 
                            
            }); 

            piApproximation = 4 * ((double)inCircle / (double)iterations);

            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));
            
        }

        private static void MonteCarloPiApproximationSerialSimulation()
        {
            double piApproximation = 0;
            int total = 0;
            int inCircle = 0; 
            double x,y = 0;
            Random rnd = new Random(); 

            while (total < iterations)
            {
                x = rnd.NextDouble(); 
                y = rnd.NextDouble();

                if ((Math.Sqrt(x*x+y*y) <= 1.0))
                    inCircle++;

                total++;                
                piApproximation =  4 * ((double)inCircle / (double)total); 
            } //while 


            Console.WriteLine();
            Console.WriteLine("Approximated Pi = {0}", piApproximation.ToString("F8"));

        }




    }
}

Data block types in Task Parallel Library

The following code sample shows data block types in Task Parallel Libary (TPL). The code is written in C# and can be run in a simple Console Application. Dataflow in TPL makes it possible to build actor-based programming and orchestrate coarse-grained dataflow and pipelining tasks, maintaining robustness and supporting concurrency-enabled applications. This makes it easier to construct high-performance, low latency systems.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using System.Threading.Tasks.Dataflow;

namespace TplDataFlowSimpleTests
{
    class Program
    {

        static void Main(string[] args)
        {
            ActionBlockPostingAndFaulting();

            //BufferBlockPostingAndReceiving();

            //BroadCastPostAndMultipleReceive();

            //WriteOnceBlockParallellPostsAndSingleReceive();

            //ActionBlockPostingCompleting();

            //TransformBlockPostingAndProducingOutput(); 

            //TransformManyBlockPostingsAndReceive(); 

            //BatchBlockSeveralBatches();

            //JoinBlockAritmeticCombinations(); 

            //BatchedJoinBlockPropagatingValuesAndException();

            Console.WriteLine("Press the any key to continue ..");
            Console.ReadKey();

        }

        private static void BatchedJoinBlockPropagatingValuesAndException()
        {
            Func<int, int> doWork = n =>
            {
                if (n < 0)
                    throw new ArgumentOutOfRangeException();
                return n;
            };

            var batchedJoinBlock = new BatchedJoinBlock<int, Exception>(7);

            foreach (int i in new int[] {5, 6, -7, -22, 13, 55, 0})
            {
                try
                {
                    batchedJoinBlock.Target1.Post(doWork(i));
                }
                catch (ArgumentOutOfRangeException e)
                {
                    batchedJoinBlock.Target2.Post(e);
                }
            }

            var results = batchedJoinBlock.Receive();

            foreach (int n in results.Item1)
            {
                Console.WriteLine(n);
            }

            foreach (Exception e in results.Item2)
            {
                Console.WriteLine(e.Message);
            }
        }

        private static void JoinBlockAritmeticCombinations()
        {
            var joinBlock = new JoinBlock<int, int, char>(new GroupingDataflowBlockOptions{ Greedy = true});

            joinBlock.Target1.Post(3);
            joinBlock.Target1.Post(6);

            joinBlock.Target2.Post(5);
            joinBlock.Target2.Post(4);

            joinBlock.Target3.Post('+');
            joinBlock.Target3.Post('-');

            for (int i = 0; i < 2; i++)
            {
                var data = joinBlock.Receive();

                switch (data.Item3)
                {
                    case '+':
                        Console.WriteLine("{0} + {1} = {2}", data.Item1, data.Item2, data.Item1 + data.Item2);
                        break;
                    case '-':
                        Console.WriteLine("{0} - {1} = {2}", data.Item1, data.Item2, data.Item1 - data.Item2);
                        break;
                    default:
                        Console.WriteLine("Unknown operator '{0}'.", data.Item3);
                        break;
                } //switch 
            } //for 

        }

        private static void BatchBlockSeveralBatches()
        {
            var batchBlock = new BatchBlock<int>(10);

            for (int i = 0; i < 13; i++)
            {
                batchBlock.Post(i);
            }

            batchBlock.Complete();


            Console.WriteLine("The elements of the first batch are: [{0}] ", string.Join(",", batchBlock.Receive()));
            Console.WriteLine("The elements of the second batch are: [{0}]", string.Join(",", batchBlock.Receive()));
        }

        private static void TransformManyBlockPostingsAndReceive()
        {
            var transformManyBlock = new TransformManyBlock<string, char>(s => s.ToCharArray());

            transformManyBlock.Post("Hello");
            transformManyBlock.Post("World");

            for (int i = 0; i < ("Hello" + "World").Length; i++)
            {
                Console.WriteLine(transformManyBlock.Receive());
            }
        }

        private static void TransformBlockPostingAndProducingOutput()
        {
            var transformBlock = new TransformBlock<int, double>(n => Math.Sqrt(n));

            transformBlock.Post(10);
            transformBlock.Post(20);
            transformBlock.Post(30);

            for (int i = 0; i < 3; i++)
            {
                Console.WriteLine(transformBlock.Receive());
            }
        }

        private static void ActionBlockPostingCompleting()
        {
            var cts = new CancellationTokenSource();
            var actionBlock = new ActionBlock<int>(n =>
            {
                Console.WriteLine(n);
                Thread.Sleep(1000);
                if (n > 50)
                    cts.Cancel();
            }, new ExecutionDataflowBlockOptions{ MaxDegreeOfParallelism = 2, MaxMessagesPerTask = 1, CancellationToken = cts.Token });

            for (int i = 0; i < 10; i++)
            {
                actionBlock.Post(i*10); 
         
            }

            actionBlock.Complete();
            try
            {
                actionBlock.Completion.Wait(cts.Token);
            }
            catch (OperationCanceledException ex)
            {
                Console.WriteLine(ex.Message);
            }
        }

        private static void WriteOnceBlockParallellPostsAndSingleReceive()
        {
            var writeOnceBlock = new WriteOnceBlock<string>(null);

            Parallel.Invoke(
                () => writeOnceBlock.Post("Message 1"),
                () => writeOnceBlock.Post("Message 2"),
                () => writeOnceBlock.Post("Message 3"));

            Console.WriteLine(writeOnceBlock.Receive());
        }

        private static void BroadCastPostAndMultipleReceive()
        {
            var broadcastBlock = new BroadcastBlock<double>(null);

            broadcastBlock.Post(Math.PI);

            for (int i = 0; i < 3; i++)
            {
                Console.WriteLine(broadcastBlock.Receive());
            }
        }

        private static void BufferBlockPostingAndReceiving()
        {
            var bufferBlock = new BufferBlock<int>();

            for (int i = 0; i < 3; i++)
            {
                bufferBlock.Post(i);
            }

            for (int i = 0; i < 4; i++)
            {
                try
                {
                    Console.WriteLine(bufferBlock.Receive(new TimeSpan(0, 0, 2)));
                }
                catch (TimeoutException tie)
                {
                    Console.WriteLine("Exception of type: {0} with message: {1}", tie.GetType().Name, tie.Message);
                }
            }
        }

        private static void ActionBlockPostingAndFaulting()
        {
            var throwIfNegative = new ActionBlock<int>(n =>
            {
                Console.WriteLine("n = {0}", n);
                if (n < 0)
                    throw new ArgumentOutOfRangeException();
            });

            throwIfNegative.Completion.ContinueWith(
                task => { Console.WriteLine("The status of the completion task is '{0}'", task.Status); });

            throwIfNegative.Post(0);
            throwIfNegative.Post(-1);
            throwIfNegative.Post(1);
            throwIfNegative.Post(2);

            throwIfNegative.Complete();

            try
            {
                throwIfNegative.Completion.Wait();
            }
            catch (AggregateException ae)
            {
                ae.Handle(e =>
                {
                    Console.WriteLine("Encountered {0}: {1}", e.GetType().Name, e.Message);
                    return true;
                });
            }

          
        }
    }
}

Producer-consumer scenario with BlockingCollection of Task Parallell Library

The BlockingCollection of Task Parallel Library or TPL supports Producer-Consumer scenarios well. This is a scenario or pattern, where you want to start processing items as soon as they are available. BlockingCollection makes this easy, as long as you follow the convention. In this article, a simple scenario with five producers and one consumer is presented. We create an array of tasks (Task[]) that we wait on, using the Task.WaitAll(Task[]) method. This creates a barrier in our main method of which we call the CompleteAdding() method of BlockingCollection to signal that we have no more items to add. The consumer, which here is a single Task, will use standard foreach iteration and call the method GetConsumingEnumerable() on the BlockingCollection to get the collection to iterate over. Note that we will be inside the foreach loop until the CompleteAdding() method is called on the BlockingCollection, i.e. the iteration is halted until CompleteAdding() is called and no more items are available.

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace TestBlockingCollectionSecond
{

    public class Node
    {

        public int ManagedThreadId { get; set; }

        public int Value { get; set; }

        public override string ToString()
        {
            return string.Format("Inside ManagedThreadId {0}, Value: {1}", ManagedThreadId, Value);
        }

    }

    class Program
    {
        static void Main(string[] args)
        {

            BlockingCollection<Node> bc = new BlockingCollection<Node>();

            Console.WriteLine("Producer-consumer scenario BlockingCollection.\nFive producers, one consumer scenario");

            var rnd = new Random();

            var producers = new List<Task>();

            for (int i = 0; i < 5; i++)
            {

                var producer = Task.Factory.StartNew(() =>
                {                  
                    for (int j = 0; j < 5; j++)
                    {
                        Thread.Sleep(rnd.Next(100,300));
                        bc.Add(new Node { ManagedThreadId = Thread.CurrentThread.ManagedThreadId, Value = rnd.Next(1, 100) });
                    }
                });

                producers.Add(producer);

            }

            var consumer = Task.Factory.StartNew(() =>
            {
                foreach (Node item in bc.GetConsumingEnumerable())
                {
                    Console.WriteLine(item);
                }

                Console.WriteLine("Consumer all done!");             
            });

            Task.WaitAll(producers.ToArray());

            bc.CompleteAdding();
       

            Console.WriteLine("Hit any key to exit program");
            Console.ReadKey();         


        }
    }
}

If you are not sure when to signal CompleteAdding(), you can keep taking items from the BlockingCollection, using one single consumer or multiple, but remember to catch InvalidOperationException, in case there are no more items available, that is - after CompleteAdding() method has been called on the BlockingCollection. Note that the while loop below is wrapped in a consumer Task that is once again started before the Task.WaitAll call on the producers array, to start consuming right away. Our exit condition here is the flag bc.IsAddingCompleted is set to true, i.e. a call to CompleteAdding is performed. The benefit of using the GetConsumingEnumerable() here is having not to deal with exceptions and boolean flag of completed adding items in the consumer Block, since calling the Take() method on the collection after the CompleteAdding() method is called will throw an InvalidOperationException.

            var consumer = Task.Factory.StartNew(() => {

                try{
                while (!bc.IsAddingCompleted){
                    Console.WriteLine("BlockingCollection element: " + bc.Take());
                }
                }
                Console.WriteLine("Consumer all done!");
                catch (InvalidOperationException ioe){
                    //Console.WriteLine(ioe.Message);
                }
            });