/*

	Program to calculate ruin-theory approach to mortality.

	See "Applying survival models to pensioner mortality data",	details of
	which can be found at http://www.richardsconsulting.co.uk/laws.html

	Copyright (c) 2008 Stephen Richards.  All rights reserved.

*/

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

#define COMMAND 0
#define LAMBDA  1
#define MU      2
#define R       3
#define X       4
#define N       5
#define ARGC    6

#define AGE_LIMIT 120	/*	Maximum age.	*/

int gidum = -1,						// For random-number generators.
    giSeed = -1;					// Copy of initial seed for random-number generators.

/*	Function to return minimum of two real values.	*/
double minimum(const double a, const double b)
{
  if (a < b)
    return a;

  return b;
}

// Park and Miller's random-number generator with Bays-Durham shuffle
// Implementation inspired by Press et al (2002)
double U(void)
{
  const int IA = 16807, IM = 2147483647, IQ = 127773, IR = 2836, MASK = 123459876;
  const double AM = 1.0 / (double)IM;
  int k;
  double ans;

  gidum ^= MASK;
  k = gidum / IQ;
  gidum = IA * (gidum - k * IQ) - IR * k;
  if (gidum < 0)
    gidum += IM;

  ans = AM * gidum;
  gidum ^= MASK;
  return ans;
}

int main(const int argc, const char * argv[])
{
  int rc = EXIT_FAILURE;

  if (argc != ARGC)
    fprintf(stderr, "Usage: ruin lambda mu r x n\n");
  else
  {
    const double lambda = atof(argv[LAMBDA]),
                 mu     = atof(argv[MU    ]),
                 r      = atof(argv[R     ]),
                 x      = atof(argv[X     ]);
    const unsigned long int n = atol(argv[N]);

    if (lambda <= 0.0)
      fprintf(stderr, "lambda must be greater than zero!\n");
    else if (mu <= 0.0)
      fprintf(stderr, "mu must be greater than zero!\n");
    else if (r < 0.0)
      fprintf(stderr, "r must be non-negative!\n");
    else if (x <= 0.0)
     fprintf(stderr, "x must be greater than zero!\n");
    else if (n <= 0)
     fprintf(stderr, "n must be greater than zero!\n");
    else
    {
      double dLifetime, dLifeLived[AGE_LIMIT+1], ux, t, damage;
      unsigned long int i, j, Deaths[AGE_LIMIT+1], Lives[AGE_LIMIT+1], ageatend;
      int Dead;

      /*	Initialise counters.	*/
      for (i=0; i<=AGE_LIMIT; i++)
      {
        Deaths[i] =
        Lives [i] = 0;
        dLifeLived[i] = 0.0;
      }

      for (i=0; i<n; i++)
      {
        dLifetime = 0.0;
        ux = x;
        Dead = 0;
        while ( (Dead==0) && (dLifetime < AGE_LIMIT) )
        {
          t = -log(U()) / lambda;		/*	Simulate time to next event.	*/
          ux = minimum(ux + t * r, x);	/*	Recover to maximum of x.		*/
          dLifetime += t;
          damage = -log(U()) / mu;			/*	Simulate damage.				*/
          ux -= damage;
          if (ux < 0.0)
            Dead = 1;
		}

        ageatend = (unsigned long int)floor(dLifetime);
        if ( (Dead == 1) && (ageatend <= AGE_LIMIT) )
          Deaths[ageatend]++;

        for (j=0; j<=minimum(ageatend, AGE_LIMIT); j++)		/*	Increment life counts for all ages started.			*/
          Lives[j]++;

        for (j=0; j<minimum(ageatend, AGE_LIMIT); j++)		/*	Add whole years of life for all ages except last.	*/
          dLifeLived[j] += 1.0;

        if (ageatend <= AGE_LIMIT)		/*	Add any last whole or fractional year of life.		*/
          dLifeLived[ageatend] += dLifetime - ageatend;

        if ((i % 100) == 0)
          fprintf(stderr, "%8ld\b\b\b\b\b\b\b\b", n - i);
      }

      /*	Output crude estimates of force of mortality at age i+0.5 */
      printf("lambda=%f, mu=%f, r=%f, x=%f, n=%d\n", lambda, mu, r, x, n);
      printf("  Age    Mu    Deaths   Lives Exposure\n");
      for (i=0; i<=AGE_LIMIT; i++)
        if (dLifeLived[i] > 0.0)
          printf("%5.1f %7.5f %7d %7d %f\n", i+0.5, Deaths[i] / dLifeLived[i], Deaths[i], Lives[i], dLifeLived[i]);

      rc = EXIT_SUCCESS;
    }
  }

  return rc;
}