General Shock Model Program

I attached source code for a generalized command-line driven program that takes input discovery data (+ other parameters) and which then generates production curves from the date of first discovery to the present time (and on out into the future). The code basically refactors the specific code from other posts in the Oil Shock Model series, to get a better separation of data from the algorithms.

The program reads the input parameters from either a set of files (for yearly data) or from the term following the option (for fixed constants). It doesn't matter what units you use as long as they remain consistent with the other ones. The program sends the output data to the shell terminal.

-d

File containing discovery data (continuous years, one entry per line)

-s

File containing shock data (only years with discontinuities), two entries per line, "year" + (space) + "extraction rate that year". The program interpolates between entries, so to prevent shocks, give it only two lines (start and finish times) with the same extraction rate on each line.

-p

Optional file containing production data (continuous years, one entry per line). The program doesn't do anything with this data; it only carries it over to make it convenient for comparison to the estimated production.

-1

The 1/e rate for remaining fallow (i.e. 1 over the average time between discovery and construction)

-2

The 1/e rate for finishing construction (i.e. 1 over the build time)

-3

The 1/e rate for reaching maturing (i.e. 1 over the startup time)

-t

Timespan for complete simulation

Typical invocation:

shockmodel.exe -1 0.2 -2 0.2 -3 0.2 -t 600 -d strikes.dat -p prod.dat -s shocks.dat > output

with Text_IO;
with Ada.Numerics.Elementary_Functions;
with GNAT.Command_Line;
with Ada.Float_Text_IO;

procedure ShockModel is
  type Flt is array (Natural range <>) of Float;
  type Flt_Access is access all Flt;

  Expansion_Factor : constant := 10; -- Year subtics

  Time_Span : Integer := 600;

  -- Safe array retrieval function
  function Get (Arr : Flt; I : Integer) return Float is
  begin
    if I in Arr'Range then
      return Arr (I);
    else
      return 0.0;
    end if;
  end Get;

  -- Discretized convolution function
  function "*" (A, B : in Flt) return Flt is
    Total : constant Natural := A'Length + B'Length;
    C     : Flt (0 .. Total);
    V     : Float;
  begin
    for J in  0 .. Total loop
      V := 0.0;
      for I in  0 .. J loop
        V := V + Get (A, I) * Get (B, J - I);
      end loop;
      C (J) := V;
    end loop;
    return C;
  end "*";

  -- Discretized exponential function
  function Exponential (L : Natural; Alpha : Float) return Flt is
    use Ada.Numerics.Elementary_Functions;
    R : Flt (0 .. L - 1);
  begin
    for I in  0 .. L - 1 loop
      R (I) := Alpha * Exp (-Alpha * Float (I));
    end loop;
    return R;
  end Exponential;

  -- Fills in data between discrete years
  function Expander (Arr : Flt) return Flt is
    R : Flt (0 .. Arr'Length * Expansion_Factor - 1);
  begin
    for I in  R'Range loop
      R (I) := Arr (I / Expansion_Factor);
    end loop;
    return R;
  end Expander;

  -- Safe function for reading forcing function
  function Discovery_Window
   (Discovery : Flt;
    Year      : Float;
    Start     : Float)
    return      Float
  is
    Y : constant Float := Year - Start;
  begin
    if Y < 0.0 then
      return 0.0;
    end if;
    return Discovery (Integer (Y) * Expansion_Factor);  -- Subdivide by factor
  end Discovery_Window;

  type Shock is record
    Year : Float; -- Year of shock
    Rate : Float; -- Rate of extraction
  end record;
  type Shocks is array (Natural range <>) of Shock;
  type Shocks_Access is access all Shocks;

  -- Make the shock a continuous function
  -- by interpolating between points in the list
  function Interpolate_Shocks (Year : Float; S : Shocks) return Float is
    K : Integer := S'Last - 1;
  begin
    for J in  S'First + 1 .. S'Last loop
      if S (J).Year > Year then
        K := J - 1;
        exit;
      end if;
    end loop;
    return S (K).Rate +
           (S (K + 1).Rate - S (K).Rate) * (Year - S (K).Year) /
           (S (K + 1).Year - S (K).Year);
  end Interpolate_Shocks;

  -- Strikes into Mature Discovery
  Strikes   : Flt_Access;
  Discovery : Flt_Access;

  Phases : Flt (1 .. 3) := (others => 0.0);

  -- Generates the mature discovery by repeated convolution
  procedure Phasing is
    Avg_Fallow : constant Flt :=
      Exponential (Time_Span, Phases (1) / Float (Expansion_Factor));
    Avg_Build  : constant Flt :=
      Exponential (Time_Span, Phases (2) / Float (Expansion_Factor));
    Avg_Mature : constant Flt :=
      Exponential (Time_Span, Phases (3) / Float (Expansion_Factor));
  begin
    Discovery :=
     new Flt'(Expander (Strikes.all) * Avg_Fallow * Avg_Build * Avg_Mature);
  end Phasing;

  S          : Shocks_Access;
  Production : Flt_Access;  -- To compare against, optional

  procedure Compute is
    C      : Float          := 0.0;  -- Reserve
    Rate   : Float;         -- Extraction Rate
    P      : Float;         -- Production Rate
    DT     : constant Float := 0.01; -- 1/100th of year
    Start  : constant Float := S (S'First).Year;
    Finish : constant Float := S (S'Last).Year;
    Time   : Float          := Start;
    V      : Float          := 0.0; -- Volume of oil extracted
    Data   : Float          := 0.0;
    PIndex : Integer;
  begin
    loop
    -- Integration of discovery window against extraction rate
      Rate := Interpolate_Shocks (Time, S.all);
      C    := C +
              Discovery_Window (Discovery.all, Time, Start) * DT -
              C * Rate * DT;
      P    := Rate * C;
      V    := V + Rate * C * DT;
      if Time > Float (Integer (Time)) - DT / 2.0 and
         Time < Float (Integer (Time)) + DT / 2.0
      then  -- Print output only once per year
        PIndex := Integer (Time - Start);
        if PIndex in Production'Range then
          Data := Production (PIndex);
        else
          Data := 0.0;
        end if;
        Text_IO.Put_Line
         (Integer'Image (Integer (Time)) & "," & P'Img & "," & Data'Img);
      end if;
      Time := Time + DT;
      exit when Time > Finish;
    end loop;
  end Compute;

  function Get_Flt_Data (File : in String) return Flt_Access is
    FT  : Text_IO.File_Type;
    function Get_Data (Data : in Flt) return Flt is
      Value : Flt (0 .. 0);
    begin
      Ada.Float_Text_IO.Get (FT, Value (0));
      Text_IO.Skip_Line (FT);
      return Get_Data (Data & Value);
    exception
      when Text_IO.End_Error =>
        return Data;
    end Get_Data;
  begin
    Text_IO.Open (FT, Text_IO.In_File, File);
    return new Flt'(Get_Data (Flt'(1 .. 0 => 0.0)));
  end Get_Flt_Data;

  function Get_Shock_Data (File : in String) return Shocks_Access is
    FT  : Text_IO.File_Type;
    function Get_Data (Data : in Shocks) return Shocks is
      Pair : Shocks (0 .. 0);
    begin
      Ada.Float_Text_IO.Get (FT, Pair (0).Year);
      Ada.Float_Text_IO.Get (FT, Pair (0).Rate);
      Text_IO.Skip_Line (FT);
      return Get_Data (Data & Pair);
    exception
      when Text_IO.End_Error =>
        return Data;
    end Get_Data;
  begin
    Text_IO.Open (FT, Text_IO.In_File, File);
    return new Shocks'(Get_Data (Shocks'(1 .. 0 => (0.0, 0.0))));
  end Get_Shock_Data;

begin
  -- Command line options
  loop
    case GNAT.Command_Line.Getopt ("1: 2: 3: t: d: s: p:") is
    when ASCII.NUL =>
      exit;
    when '1' =>
      Phases (1) := Float'Value (GNAT.Command_Line.Parameter);
    when '2' =>
      Phases (2) := Float'Value (GNAT.Command_Line.Parameter);
    when '3' =>
      Phases (3) := Float'Value (GNAT.Command_Line.Parameter);
    when 't' =>
      Time_Span := Integer'Value (GNAT.Command_Line.Parameter);
    when 'd' =>
      Strikes := Get_Flt_Data (GNAT.Command_Line.Parameter);
    when 's' =>
      S := Get_Shock_Data (GNAT.Command_Line.Parameter);
    when 'p' =>
      Production := Get_Flt_Data (GNAT.Command_Line.Parameter);
    when others =>
      raise Program_Error;         -- cannot occur!
    end case;
  end loop;
  Phasing;
  Compute;
exception
  when GNAT.Command_Line.Invalid_Switch =>
    Text_IO.Put_Line ("Invalid Switch " & GNAT.Command_Line.Full_Switch);
  when GNAT.Command_Line.Invalid_Parameter =>
    Text_IO.Put_Line ("No parameter for " & GNAT.Command_Line.Full_Switch);
end ShockModel;