Re: Benchmark Ada, please

[Guillaume Foliard <guifo@wanadoo.fr>]

Victor Porton wrote:

> Somebody, write an Ada benchmark for this comparison of programming
> languages:
> 
> https://github.com/jesusfv/Comparison-Programming-Languages-Economics

Here is it:

---------------------------------------------------------------------
with Ada.Execution_Time;
with Ada.Numerics.Long_Elementary_Functions;
with Ada.Real_Time;
with Ada.Text_IO;

use Ada.Numerics.Long_Elementary_Functions;

use type Ada.Execution_Time.CPU_Time;

procedure Rbc_Ada
is
   
   Grid_Capital_Length      : constant := 17820;
   Grid_Productivity_Length : constant := 5;
   
   type Grid_Capital_Index is 
      new Integer range 0 .. Grid_Capital_Length - 1;
   
   type Grid_Productivity_Index is
      new Integer range 0 .. Grid_Productivity_Length - 1;
   
   type Grid_Array_Type is
      array (Grid_Capital_Index, Grid_Productivity_Index)
      of Long_Float;
   
   Null_Grid : constant Grid_Array_Type := (others => (others => 0.0));
   
   -- 1. Calibration
   
   -- Elasticity of output w.r.t. capital
   
   Alpha : constant Long_Float := 0.33333333333;
   
   -- Discount factor
   
   Beta : constant Long_Float := 0.95;
   
   -- Productivity values
   
   Productivities : constant
      array (Grid_Productivity_Index) 
      of Long_Float :=
      (0.9792, 0.9896, 1.0000, 1.0106, 1.0212);
   
   -- Transition matrix
   
   Transition : constant
      array (Grid_Productivity_Index, Grid_Productivity_Index) 
      of Long_Float :=
      ((0.9727, 0.0273, 0.0000, 0.0000, 0.0000),
       (0.0041, 0.9806, 0.0153, 0.0000, 0.0000),
       (0.0000, 0.0082, 0.9837, 0.0082, 0.0000),
       (0.0000, 0.0000, 0.0153, 0.9806, 0.0041),
       (0.0000, 0.0000, 0.0000, 0.0273, 0.9727));
   
   -- 2. Steady State
   
   Capital_Steady_State : constant Long_Float :=
      (Alpha * Beta) ** (1.0 / (1.0 - Alpha));
   
   Output_Steady_State : constant Long_Float := 
      Capital_Steady_State ** Alpha;
   
   Consumption_Steady_State : constant Long_Float := 
      Output_Steady_State - Capital_Steady_State;

   Grid_Capital_Next_Index : Grid_Capital_Index;
   
   Grid_Capital : array (Grid_Capital_Index) of Long_Float := 
      (others => 0.0);
  
   Output                  : Grid_Array_Type := Null_Grid;
   Value_Function          : Grid_Array_Type := Null_Grid;
   Value_Function_New      : Grid_Array_Type := Null_Grid;
   Policy_Function         : Grid_Array_Type := Null_Grid;
   Expected_Value_Function : Grid_Array_Type := Null_Grid;

   Max_Difference   : Long_Float := 10.0;
   Diff             : Long_Float;
   Diff_High_So_Far : Long_Float;
   
   Tolerance : constant := 0.0000001;
   
   Value_High_So_Far : Long_Float;
   Value_Provisional : Long_Float;
   Consumption       : Long_Float;
   Capital_Choice    : Long_Float;
   Iteration         : Integer := 0;
   Cpu_Time_Start    : Ada.Execution_Time.CPU_Time;
   Cpu_Time_End      : Ada.Execution_Time.CPU_Time;
begin
   Cpu_Time_Start := Ada.Execution_Time.Clock;
   
   Ada.Text_IO.Put_Line
      ("Output =" & Output_Steady_State'Img
       & ", Capital =" & Capital_Steady_State'Img
       & ", Consumption =" & Consumption_Steady_State'Img);
	
   -- We generate the grid of capital

   for Index in Grid_Capital'Range
   loop
      Grid_Capital (Index) :=
         0.5 * Capital_Steady_State + 0.00001 * Long_Float (Index);
   end loop;

   -- We pre-build output for each point in the grid
   
   for Productivity_Index in Grid_Productivity_Index
   loop
      for Capital_Index in Grid_Capital_Index
      loop
         Output (Capital_Index, Productivity_Index) := 
            Productivities (Productivity_Index)
            * Grid_Capital (Capital_Index) ** Alpha;
      end loop;
   end loop;

   -- Main iteration

   while Max_Difference > Tolerance
   loop
      for Productivity_Index in Grid_Productivity_Index   
      loop
         for Capital_Index in Grid_Capital_Index  
         loop
            Expected_Value_Function (Capital_Index, Productivity_Index) :=
               0.0;
            
            for Productivity_Next_Index in Grid_Productivity_Index   
            loop
               Expected_Value_Function (Capital_Index, Productivity_Index) 
:=
                  Expected_Value_Function (Capital_Index, 
Productivity_Index)
                  + Transition (Productivity_Index, Productivity_Next_Index)
                  * Value_Function (Capital_Index, Productivity_Next_Index);
            end loop;
         end loop;
      end loop;
      
      for Productivity_Index in Grid_Productivity_Index   
      loop
         -- We start from previous choice (monotonicity of policy function)
         
         Grid_Capital_Next_Index := 0;
         
         for Capital_Index in Grid_Capital_Index
         loop
            Value_High_So_Far := -100000.0;
            Capital_Choice    := Grid_Capital (0);
            
            for Capital_Next_Index in 
               Grid_Capital_Next_Index .. Grid_Capital_Index'Last
            loop
               Consumption := 
                  Output (Capital_Index, Productivity_Index)
                  - Grid_Capital (Capital_Next_Index);
               
               Value_Provisional :=
                  (1.0 - Beta) * Log (Consumption)
                  + Beta * Expected_Value_Function (Capital_Next_Index,
                                                    Productivity_Index);
               
               if Value_Provisional > Value_High_So_Far
               then
                  Value_High_So_Far := Value_Provisional;
                  Capital_Choice := Grid_Capital (Capital_Next_Index);
                  Grid_Capital_Next_Index := Capital_Next_Index;
                  
               else
                  exit;
               end if;
               
               Value_Function_New (Capital_Index, Productivity_Index) := 
                  Value_High_So_Far;
               
               Policy_Function (Capital_Index, Productivity_Index) :=
                  Capital_Choice;
            end loop;
         end loop;
      end loop;
      
      Diff_High_So_Far := -100000.0;
      
      for Productivity_Index in Grid_Productivity_Index   
      loop
         for Capital_Index in Grid_Capital_Index
         loop
            Diff := 
               abs (Value_Function (Capital_Index, Productivity_Index)
                   - Value_Function_New (Capital_Index, 
                                        Productivity_Index));
 
            if Diff > Diff_High_So_Far
            then
               Diff_High_So_Far := Diff;
            end if;
            
            Value_Function (Capital_Index, Productivity_Index)
               := Value_Function_New (Capital_Index, Productivity_Index);
         end loop;
      end loop;
      
      Max_Difference := Diff_High_So_Far;

      Iteration := Iteration + 1;
      
      if Iteration mod 10 = 0 or Iteration = 1
      then
         Ada.Text_IO.Put_Line ("Iteration =" & Iteration'Img
                               & ", Sup Diff =" & Max_Difference'Img);
      end if;
   end loop;
   
   Ada.Text_IO.Put_Line ("Iteration =" & Iteration'Img
                         & ", Sup Diff =" & Max_Difference'Img);
   Ada.Text_IO.New_Line;
   Ada.Text_IO.Put_Line ("My check =" & Policy_Function (999, 2)'Img);
   Ada.Text_IO.New_Line;   

   Cpu_Time_End := Ada.Execution_Time.Clock;
   
   Ada.Text_IO.Put_Line
      ("Elapsed time is ="
       & Ada.Real_Time.To_Duration (Cpu_Time_End - Cpu_Time_Start)'Img);
end Rbc_Ada;
---------------------------------------------------------------------

This is mostly a line to line translation from RBC_CPP.cpp. I have
added a few type declarations though.

> It seems that C++ was the fastest (faster than Fortran), but Ada may be
> even faster.

Here are the numbers with GNAT GPL 2014 on a Core2 Q9650 @ 3.00GHz:

$ gnatmake -O3 rbc_ada.adb
$ time ./rbc_ada
...
Elapsed time is = 1.966112682


As for the C++ version:

$ g++ -o testc -O3 RBC_CPP.cpp
$ time ./testc
... 
Elapsed time is   = 3.12033

So the Ada version is significantly faster. I suppose it is mainly because 
the Ada compiler has vectorized more loops than the C++ compiler (add -
ftree-vectorizer-verbose=2 to the above compilation commands to check by 
yourself).

> If we succeed, we would advertise Ada as the fastest(!) programming
> language (after assembler).

Feel free to advertise, using this Ada code as you wish.

-- 
Guillaume Foliard