/* 
 * Template of actuarial calculation of future cash flow for a group of 
 * life policies. Group of policies are described by:
 *
 *  number_of_policies
 *  length
 *  premium
 *  beneficiary_return_factor
 *  initial_age                  
 *
 *  This is to mimic Alef software behaviour and try different ways of
 *  improving performances.
 *
 * Author: M. Cestari
 * Year: 2014  
*/

#define MAN 0
#define WOMAN 1

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_rng.h> 
#include <string.h>    
#include "cclock.h"

struct PolicyParams {
    int length; /* length in years*/  
    int premium; /* monthly premium in Euros */   
    int beneficiary_return_factor; /* to calculate the return value for the beneficiary */ 
    int initial_age; /* initial age of the insured */                
};

int read_input(char *name, int *number, struct PolicyParams *params) {

    /* Reading settings of the group of policies.  
     * Store the results on a PolicyParams structure */ 

    FILE *in_file;  /* file containing the input data */
    char line[100]; /* line of 100 characters*/  
    int scan_count; /* Number of parameters scanned */

    in_file = fopen(name, "r");
    if (in_file == NULL) { /* Test for error */
       fprintf(stderr, "Error: Unable to input file %s\n",name);
       exit (8);
    }
    
    /* read number of policies*/
    fgets(line, sizeof(line), in_file);
    scan_count = sscanf(line,"%d", number);

    /* read policy length*/
    fgets(line, sizeof(line), in_file);
    scan_count = sscanf(line,"%d", &(params->length));

    /* read premium*/
    fgets(line, sizeof(line), in_file);
    scan_count = sscanf(line,"%d", &(params->premium));
    
    /* read beneficiary_return_factor*/
    fgets(line, sizeof(line), in_file);
    scan_count = sscanf(line,"%d", &(params->beneficiary_return_factor));
    
    /* read initial_age*/
    fgets(line, sizeof(line), in_file);
    scan_count = sscanf(line,"%d", &(params->initial_age));

    fclose(in_file);
}

float death_probability(int updated_age, int gender){
    float actuarial_life_table_man[] = { 
	    0.006990, 0.000447, 0.000301, 0.000233, 0.000177,
	    0.000161, 0.000150, 0.000139, 0.000123, 0.000105,
	    0.000091, 0.000096, 0.000135, 0.000217, 0.000332,
	    0.000456, 0.000579, 0.000709, 0.000843, 0.000977,
	    0.001118, 0.001250, 0.001342, 0.001382, 0.001382, 
	    0.001370, 0.001364, 0.001362, 0.001373, 0.001393,
	    0.001419, 0.001445, 0.001478, 0.001519, 0.001569,
	    0.001631, 0.001709, 0.001807, 0.001927, 0.002070,
            0.002234, 0.002420, 0.002628, 0.002860, 0.003117,  
            0.003396, 0.003703, 0.004051, 0.004444, 0.004878, 
            0.005347, 0.005838, 0.006337, 0.006837, 0.007347,
            0.007905, 0.008508, 0.009116, 0.009723, 0.010354,
            0.011046, 0.011835, 0.012728, 0.013743, 0.014885,
            0.016182, 0.017612, 0.019138, 0.020752, 0.022497,
            0.024488, 0.026747, 0.029212, 0.031885, 0.034832,
            0.038217, 0.042059, 0.046261, 0.050826, 0.055865,
            0.061620, 0.068153, 0.075349, 0.083230, 0.091933,
            0.101625, 0.112448, 0.124502, 0.137837, 0.152458,
            0.168352, 0.185486, 0.203817, 0.223298, 0.243867, 
            0.264277, 0.284168, 0.303164, 0.320876, 0.336919,
            0.353765, 0.371454, 0.390026, 0.409528, 0.430004,
            0.451504, 0.474079, 0.497783, 0.522673, 0.548806,
            0.576246, 0.605059, 0.635312, 0.667077, 0.700431,
            0.735453, 0.772225, 0.810837, 0.851378, 0.893947,
    };
    if (gender==MAN){
        return actuarial_life_table_man[updated_age];
    }
    else{
        return actuarial_life_table_man[updated_age];
        //return actuarial_life_table_woman[updated_age]
    }
}

int is_alive(float death_probability,  gsl_rng *r){
    /* return True if alive*/
  
    double u;

    gsl_rng_env_setup();

    u = 0.0;
    u = gsl_rng_uniform(r);
    
    if (u<=death_probability){ 
    	//printf("Random= %f  death_probability = %f \n", u, death_probability); 
        return 0;
    }
    return 1; 
}


int main(int argc, char* argv[]){

    int number_of_policies;  /* total numer of the policies in the group */
    struct PolicyParams policy_parameters; /* params of the policies group */ 

    float insured_return; 
    const int nmonths = 12; /* months per year*/
    float beneficiary_return; 
    float* cash_flow; /* array for cash flows */

    extern void test_dgemm_(int *);
    

    /* start the clock*/
    double t1 = cclock( );

    read_input(argv[1], &number_of_policies, &policy_parameters);	
 
    printf("*********************************************************************\n");
    printf("Number of considered policies: %d\n", number_of_policies);
    printf("Length of the policies: %d\n", policy_parameters.length);
    printf("Monthly Premium: %d\n", policy_parameters.premium);
    printf("Beneficiary return factor: %d\n", policy_parameters.beneficiary_return_factor);
    printf("Initial Age of the insured: %d\n", policy_parameters.initial_age);
    printf("*********************************************************************\n");


    /* we suppose an increase of insured_return value by 5% */
    insured_return = policy_parameters.premium * policy_parameters.length * nmonths;
    insured_return = insured_return * 1.05;
    beneficiary_return = insured_return * policy_parameters.beneficiary_return_factor; 

    /* calculate the number of simulation steps*/
    int number_of_steps = policy_parameters.length * nmonths;

    /* allocate array for cash_flow*/
    cash_flow = (float *)malloc((number_of_steps+1)*sizeof(float));
    for (int j=0; j<number_of_steps; j++){
        cash_flow[j] = 0.0;
    }


    int policy_completion = 0; /* True if insured outlast life policy */
    
    int premium = policy_parameters.premium;
    int initial_age = policy_parameters.initial_age;
    /* loop over the policies */
    for (int i=0; i<number_of_policies; i++){
	int gender = (2*i)/number_of_policies; // 0 Man, 1 Woman, first half Men     
        /* initialize random generator  */
	const gsl_rng_type *T; 
        gsl_rng *r;                            
        T = gsl_rng_default;
        r = gsl_rng_alloc(T);
	gsl_rng_set(r, i); // ste a different seed for each MC simulation
        /* loop over the simulation steps */
    	for (int j=0; j<number_of_steps; j++){
            int updated_age = ((initial_age*nmonths)+j)/nmonths; // every 12 months update age +1
            float death_prob = death_probability(updated_age, gender);
	    if (j%12 == 0){  // on January check whether is still alive
	        if (!is_alive(death_prob, r)){
	            //printf("age %d, step %d, death_prob %f\n", updated_age, j, death_prob);  	
                    cash_flow[j] -= beneficiary_return;
	            break; // bail out from MC loop
	        }
	    }   
	    cash_flow[j] = cash_flow[j] + premium;
	    policy_completion = (j+1)/number_of_steps; // if j+1=number_of_steps -> policy_completion
	    //int dim = 20;
	    //test_dgemm_(&dim); // this is just to waste cpu-time
	}
        gsl_rng_free(r);
	if (policy_completion) cash_flow[number_of_steps] = cash_flow[number_of_steps] - insured_return;
    }

    /* print output */
    int year;
    int trimester;
    for (int i=0; i<number_of_steps+1; i++){
	year = ((initial_age*nmonths)+i)/nmonths;
        //trimester = (4*i)/nmonths;
        printf("%3d  %3d  %15.2f\n",year, i%12+1, cash_flow[i]);  	
    }

    double t2 = cclock( );
    free(cash_flow);
    printf("*********************************************************************\n");
    printf("Elapsed time in seconds %f\n",t2-t1);   
    printf("*********************************************************************\n");

} /* end main */