/*
 * Simulation of the babysitter one-dimensional closed economy.
 * Copyright (C) 2008 Telford Tendys <telford@lnx-bsp.net>
 *
 * --------------------------------------------------------------------------------
 *
 * This program is licensed under the same terms as the Cairo library itself, which is
 * either the GNU Lesser General Public License (LGPL) version 2.1 or the Mozilla Public
 * License (MPL) version 1.1 at your option.
 * 
 * References:
 *   http://www.slate.com/id/1937/                "Baby-Sitting the Economy" by Paul Krugman
 *   http://web.mit.edu/krugman/www/babysit.html  "Baby-Sitting the Economy" by Paul Krugman
 *   http://web.mit.edu/krugman/www/howfast.html  "How Fast can the U.S. Economy Grow?" by Paul Krugman
 *   http://economistsview.typepad.com/economistsview/2006/05/babysitting_the.html
 *   http://www.galmarley.com/framesets/fs_fiat_money.htm
 *   http://krugman.blogs.nytimes.com/            Paul Krugman's "blog" at New York Times
 *
 * Notes:
 *   Closed economy because the number of couples is constant for all time
 *   (no one here gets out alive, but no one gets in either).
 *
 *   One-dimensional because there is only one commodity to offer for sale,
 *   thus no need for exchange rate pricing between commodities.
 *
 *   The system also has zero inertia and decisions are resolved instantaneously.
 *
 *   There is a guaranteed safety net for unmet demand -- people simply stay home
 *   and miss their chance to go out (no one starves for lack of babysitting).
 *
 *   Completely fixed money supply, no one accepts IOU.
 *
 * See also:
 *   http://www.cairographics.org/
 *   http://wiki.mozilla.org/Cairo
 *   http://www.gnu.org/licenses/lgpl.html
 *   http://www.mozilla.org/MPL/
 *
 * To compile, you want something like:
 *   gcc -o bs_2008-11-15 -Wall -fomit-frame-pointer -O5 bs_2008-11-15.c -lm -I /usr/include/cairo -lcairo
 */
#define _GNU_SOURCE 1
#include <cairo.h>
#include <getopt.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#define NUM_COUPLES 500
#define DEBUG_DETAILS 0
#define TOT_DAYS (365.25 * 3) /* The three year plan */

/* For the timeplot */
#define RHS_MARGIN 50
#define LHS_MARGIN 5
#define TOP_MARGIN 3
#define BOT_MARGIN 3


typedef struct state_S state_T;

typedef struct param_S
{
	double season_amplitude;
	double season_period;
	double base_hrs_mean;
	double base_hrs_sigma;
	double base_hrs_factor;
	double start_balance;
	double savings_target;
	double savings_target_factor;
	double give_up_chance;
	double rich_threshold; /* Upper limit on savings, after this people can't be bothered working */
	int base_hrs_order;
	int max_kids;
	int phonecalls; /* How many times people ring till they give up */
	int overspend_ok; /* 1 => people will deliberately go out MORE in order to spend away excess balance */
} param_T;

typedef struct couple_S
{
	double savings_target;
	double savings_target_factor;
	double balance;
	double hrs;
	int state;
	int kids_minded;
} couple_T;

struct state_S
{
	double days;
	double traded_hrs;
	double demand_hrs;
	couple_T cup[ NUM_COUPLES ];
	int (*daytots)( int *, state_T * );
};

/* -------------------------------------------------------------------------------- */

enum
{
	STATE_WANT_TO_GO_OUT          = 0,
	STATE_WANT_TO_STAY_HOME       = 1,
	STATE_INSUFFICIENT_FUNDS      = 2,
	STATE_NO_SITTER               = 3,
	STATE_GO_OUT_HAVE_FUN         = 4,
	STATE_IDLE_RICH               = 5,
	STATE_SAVING_FOR_TARGET       = 6,
};

#define NUM_STATES 7

/* -------------------------------------------------------------------------------- */

int setup_start( param_T *P, state_T *S )
{
	int i;

	S->days = 0;
	S->traded_hrs = 0;
	S->demand_hrs = 0;
	for( i = 0; i < NUM_COUPLES; i++ )
	{
		couple_T *C = S->cup + i;
/*
 * Don't start everyone EXACTLY equal because strange things happen,
 * start with a uniform distribution. Ultimately, the starting configuration
 * should not really matter but we want to allow the money supply to be large
 * enough that the majority of people turn into idle rich (just to explore that
 * dimension).
 */
		{
			double d;
			d = rand();
			d /= ( RAND_MAX / 2 );
			C->balance = d * P->start_balance;
		}

		C->savings_target = P->savings_target;
		C->savings_target_factor = P->savings_target_factor;
	}
	return( 0 );
}

int setup_random_factor( param_T *P )
{
	switch( P->base_hrs_order )
	{
		case 1: P->base_hrs_factor = P->base_hrs_sigma / 0.5773; break;
		case 2: P->base_hrs_factor = P->base_hrs_sigma / 0.8164; break;
		case 3: P->base_hrs_factor = P->base_hrs_sigma / 1.0000; break;
		case 4: P->base_hrs_factor = P->base_hrs_sigma / 1.1547; break;
		default: fprintf( stderr, "Unsupported random generator, order=%d\n", P->base_hrs_order ); exit( -1 );
	}
	return( 0 );
}

double base_random_hrs( param_T *P )
{
	int i;
	double a = 0;
	
	for( i = 0; i < P->base_hrs_order; i++ )
	{
		double x = rand();
		x -= RAND_MAX / 2;
		x /= RAND_MAX / 2;
		a += x;
	}
	a *= P->base_hrs_factor;
	return( a + P->base_hrs_mean );
}

/*
 * Chance of giving up after every failed phonecall.
 *
 * Note that each person who gives up, provides a potential sitter for
 * someone else so we do want people to give up at some stage, but we
 * don't want them to all give up at once. Ideal situation is some smart
 * arbitrator who produces a best overall result but we let it happen by
 * random chance which is probably more realistic.
 */
int give_up( param_T *P )
{
	double x = rand();
	x /= RAND_MAX;
	if( x < P->give_up_chance ) return( 1 );
	return( 0 );
}

/*
 * Seasons go as a sine wave, caused by earth rotating around the sun.
 * In summer, everyone wants to go out more, in winter they go out less.
 */
double seasonal_adj( param_T *P, state_T *S )
{
	double t = S->days;
	t /= P->season_period;
	t *= M_PI * 2;
	return( P->season_amplitude * sin( t ));
}

/*
 * If you are below your savings target then return a negative number
 * which means you will reduce the hours spent going out and thus boost
 * savings. If you are above your savings target then go out more to
 * spend the excess.
 *
 * The factor controls how seriously the couple actually takes their
 * savings target (do they even care?)
 */
double target_adj( param_T *P, couple_T *C )
{
	if( !P->overspend_ok )
		if( C->balance > C->savings_target )
			return( 0 ); /* target achieved, adjustment not required */
	return(( C->balance - C->savings_target ) * C->savings_target_factor );
}

/*
 * Pick the hours for each couple,
 * assign the states that we already are sure of such as people with no funds.
 */
int hours_round( param_T *P, state_T *S )
{
	int i;

	for( i = 0; i < NUM_COUPLES; i++ )
	{
		double natural_hrs;

		couple_T *C = S->cup + i;
		natural_hrs = base_random_hrs( P ) + seasonal_adj( P, S );
		C->hrs = natural_hrs + target_adj( P, C );
/*
 * Arguable statistic... demand hours is the "natural" demand which includes
 * factors such as season which are outside the scope of control. Add target adjustment
 * later to include "artificial" demand adjustment caused by saving/spending habits.
 */
		if( natural_hrs > 0 ) S->demand_hrs += natural_hrs;
		C->state = STATE_WANT_TO_GO_OUT;
		C->kids_minded = 0;
		if( C->hrs > 0 )
		{
/*
 * Might as well spend the last of our funds, even if it doesn't fully meet our needs.
 * This is a bit of an arguable one. Is this person satisfied or not?
 */
			if(( C->hrs > C->balance ) && ( C->balance > 0.0001 )) C->hrs = C->balance;

			if( C->hrs <= C->balance )
			{ /* Maybe will find sitter in phonecall round, keep state 0 for now */
			}
			else
			{
				C->state = STATE_INSUFFICIENT_FUNDS;
			}
		}
		else
		{
			if( C->balance > P->rich_threshold )
			{
				C->state = STATE_IDLE_RICH;
			}
			else
			{
				if( natural_hrs > 0 )
				{
					C->state = STATE_SAVING_FOR_TARGET;
				}
				else
				{
					C->state = STATE_WANT_TO_STAY_HOME;
				}
			}
		}

#if DEBUG_DETAILS
		printf( "HR[%d] => state=%d hrs=%f balance=%f\n", i, C->state, C->hrs, C->balance );
#endif
	}
	return( 0 );
}

int phonecall_round( param_T *P, state_T *S )
{
	int i, j, k;
	int offset;

	offset = rand();
	offset %= NUM_COUPLES; /* Random start point to make fair */

/*
 * For each phonecall round, visit every couple.
 * Those couples that still want to go out, make a phonecall to some
 * other couple who will offer babysitting if they are able to (why not?)
 * If the phonecall leads to a dead end, some will get discouraged and give up.
 */
	for( j = 0; j < P->phonecalls; j++ )
	{
		for( i = 0; i < NUM_COUPLES; i++ )
		{
			couple_T *C = S->cup + offset;
			couple_T *Ck;

			switch( C->state )
			{
				case STATE_WANT_TO_GO_OUT:
					for(;;)
					{
						k = rand();
						k %= NUM_COUPLES;
						if( k != offset ) break;
					}

					Ck = S->cup + k;
					switch( Ck->state )
					{
						case STATE_WANT_TO_STAY_HOME:
						case STATE_SAVING_FOR_TARGET:
						case STATE_INSUFFICIENT_FUNDS:
						case STATE_NO_SITTER:
							if( Ck->balance > P->rich_threshold )
							{ /* Can't find a sitter, but hey, not gonna help anyone else either! */
								continue;
							}

							if( Ck->kids_minded >= P->max_kids )
							{
#if DEBUG_DETAILS
								printf( "PC[%d] => reached limit of kids [%d]\n", offset, k );
#endif
								continue;
							}
							++( Ck->kids_minded );
							Ck->balance    += C->hrs;
							C->balance     -= C->hrs;
							S->traded_hrs  += C->hrs;
							C->state = STATE_GO_OUT_HAVE_FUN;

#if DEBUG_DETAILS
							printf( "PC[%d] => trade with [%d] hrs=%f\n", offset, k, C->hrs );
#endif
							goto got_it;

						default:
							if( give_up( P ))
							{
								C->state = STATE_NO_SITTER;
							}
					}
			}
got_it:;
#if DEBUG_DETAILS
			printf( "PC[%d] => state=%d hrs=%f balance=%f\n", offset, C->state, C->hrs, C->balance );
#endif

			++offset;
			if( offset >= NUM_COUPLES ) offset = 0;
		}
	}
/*
 * Anyone left in STATE_WANT_TO_GO_OUT didn't give up but we
 * ran out of rounds so in effect, they are screwed anyhow
 */
	return( 0 );
}

/* -------------------------------------------------------------------------------- */

int daytot_printer( int *ST, state_T *S )
{
	ST[ STATE_NO_SITTER ] += ST[ STATE_WANT_TO_GO_OUT ];

	printf( "%.0f\t%d\t%d\t%d\t%d\t%d\t%.3f\t%.3f\t%.4f\n", S->days,
			ST[ 1 ], ST[ 2 ], ST[ 3 ], ST[ 4 ], ST[ 5 ], 
			S->traded_hrs, S->demand_hrs, S->traded_hrs/S->demand_hrs );
	return( 0 );
}

cairo_surface_t *surface;
cairo_t *cr;
char *months[] = { "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC", "JAN", "FEB" };
char *output_file = "output.png";

int timeplot_init( param_T *P )
{
	int i;

	surface = cairo_image_surface_create( CAIRO_FORMAT_ARGB32,
										  NUM_COUPLES + LHS_MARGIN + RHS_MARGIN,
										  TOP_MARGIN + BOT_MARGIN + TOT_DAYS );
	cr = cairo_create( surface );

	cairo_set_source_rgb( cr, 0, 0, 0 );
	cairo_paint( cr );

	cairo_set_line_width( cr, 1.8 );
	cairo_set_line_cap( cr, CAIRO_LINE_CAP_ROUND );

	cairo_translate( cr, 0, TOP_MARGIN );

/*
 * Write months down the side
 */
	cairo_set_source_rgba( cr, 1, 1, 1, 1 );
	cairo_select_font_face( cr, "Sans", CAIRO_FONT_SLANT_NORMAL,
							CAIRO_FONT_WEIGHT_BOLD );

	cairo_scale( cr, 1, 3 );
	cairo_set_font_size( cr, 17 );
	for( i = 15; i < TOT_DAYS; i += 18 )
	{
		double m;
		int mm;

		m = i * 3 * 12;
		m /= P->season_period;
		mm = floor( m + 0.5 );
		mm %= 12;
		cairo_move_to( cr, 2, i );
		cairo_show_text( cr, months[ mm ]);
	}
	cairo_scale( cr, 1, 0.3333333 );

	cairo_translate( cr, RHS_MARGIN, 0 );

//	cairo_translate( cr, 320, 240 );
//	cairo_scale( cr, 300, 230 );

	return( 0 );
}

int daytot_timeplot( int *ST, state_T *S )
{
	double x1, x2;
	x1 = 0;
	x2 = ST[ 5 ];

	ST[ STATE_NO_SITTER ] += ST[ STATE_WANT_TO_GO_OUT ];

	cairo_set_source_rgba( cr, 1, 1, 0, 0.7 ); /* Got the gold, too lazy to work */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );

	x1 = x2; x2 += ST[ 1 ];
	cairo_set_source_rgba( cr, 0, 0, 1, 0.7 ); /* Got the blues, stay home, might work */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );

	x1 = x2; x2 += ST[ 6 ];
	cairo_set_source_rgba( cr, 0, 1, 1, 0.7 ); /* Staying home, saving to meet target */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );

	x1 = x2; x2 += ST[ 2 ];
	cairo_set_source_rgba( cr, 1, 0, 0, 0.7 ); /* In the Red, can't afford to go out */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );

	x1 = x2; x2 += ST[ 3 ];
	cairo_set_source_rgba( cr, 0, 0, 0, 0.7 ); /* Absence of colour, absence of babysitter */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );

	x1 = x2; x2 += ST[ 4 ];
	cairo_set_source_rgba( cr, 0, 1, 0, 0.7 ); /* Green, good to go! */
	cairo_move_to( cr, x1, S->days ); cairo_line_to( cr, x2, S->days ); cairo_stroke( cr );
	
	return( 0 );
}

int timeplot_finish( void )
{
	cairo_destroy( cr );
	cairo_surface_write_to_png( surface, output_file );
	cairo_surface_destroy( surface );
	return( 0 );
}

int collect_stats( state_T *S )
{
	int i;
	int ST[ NUM_STATES ] = {};

	for( i = 0; i < NUM_COUPLES; i++ )
	{
		couple_T *C = S->cup + i;

		ST[ C->state ] ++;
	}
	if( S->daytots )
		( *S->daytots )( ST, S );

	return( 0 );
}

/* -------------------------------------------------------------------------------- */

int standard_setup( param_T *P )
{
	P->season_amplitude       = 2;
	P->season_period          = 365.25;
	P->base_hrs_mean          = 1;
	P->base_hrs_sigma         = 4;
	P->base_hrs_order         = 3; /* rough Gaussian, short tails */
	P->start_balance          = 10;
	P->rich_threshold         = 25;
	P->savings_target         = 10;
	P->savings_target_factor  = 0.2;
	P->give_up_chance         = 0.03;
	P->max_kids               = 2;
	P->phonecalls             = 50;
	P->overspend_ok           = 0;

	setup_random_factor( P );
	return( 0 );
}

#define SWEEP_MAX 20

double sweep_step = 0.1;
double *sweep_ptr;
double *sweep_var[ SWEEP_MAX ];
double sweep_from[ SWEEP_MAX ];
double sweep_to[ SWEEP_MAX ];
int sweep_top;

/*
 * When running sweeps, handle a dimension per nesting level
 * and keep the pointers in a stack. Note that all sweeps do
 * exactly 11 samples like: [ 0, 0.1, 0.2, ... 0.9, 1.0 ]
 * Thus all sweeps are linear (actually we have an option for 101 samples
 * making really spiffy detailed graphs)
 */
int sweeper( param_T *P, int ix )
{
	double x;
	double *p;
	state_T S = {};
	int i;

	if( ix < 0 )
	{
		S.daytots = 0;
		setup_start( P, &S );
		for( i = 0; i < TOT_DAYS; i++ )
		{
			hours_round( P, &S );
			phonecall_round( P, &S );
			collect_stats( &S );
			S.days ++;
		}
		for( i = sweep_top; i--; )   /* Print one line per sweep */
		{
			p = sweep_var[ i ];
			printf( "%.3f\t", *p );
		}
		printf( "%.2f\t%.2f\t%.2f\n",
				S.demand_hrs / ( NUM_COUPLES * TOT_DAYS ),
				S.traded_hrs / ( NUM_COUPLES * TOT_DAYS ),
				100 * S.traded_hrs / S.demand_hrs );
		return( 0 );
	}

	p = sweep_var[ ix ];
	for( x = 0; x <= 1.000001; x += sweep_step )
	{
		*p = x * sweep_to[ ix ] + ( 1-x ) * sweep_from[ ix ];
		sweeper( P, ix - 1 );
	}
	return( 0 );
}

int main( int argc, char *argv[])
{
	param_T P = {};
	state_T S = {};
	int i;

	standard_setup( &P );

	for(;;)
	{
		int option_index = 0;
		int c;

		static struct option long_options[] =
			{
				{ "base-hrs-mean", 1, 0, 'b' },
				{ "base-hrs-sigma", 1, 0, 's' },
				{ "max-kids", 1, 0, 'k' },
				{ "output-file", 1, 0, 'x' },
				{ "overspend-ok", 1, 0, 'o' },
				{ "rich-threshold", 1, 0, 'r' },
				{ "savings-target", 1, 0, 't' },
				{ "savings-target-factor", 1, 0, 'f' },
				{ "season-amplitude", 1, 0, 'a' },
				{ "season-period", 1, 0, 'p' },
				{ "start-balance", 1, 0, 'm' },
				{ "sweep-fine", 0, 0, 'z' },
				{ "sweep-to", 1, 0, 'q' },
				{ 0, 0, 0, 0 }
			};

		c = getopt_long( argc, argv, "", long_options, &option_index);
		if( c < 0 ) break;
		switch( c )
		{
			case 0:
				printf( "option %s", long_options[option_index].name );
				if( optarg )
					printf( " with arg %s", optarg );
				printf( "\n" );
				break;

			case 'a': sweep_ptr = &P.season_amplitude;        goto double_var;
			case 'p': sweep_ptr = &P.season_period;           goto double_var;
			case 'b': sweep_ptr = &P.base_hrs_mean;           goto double_var;
			case 's': sweep_ptr = &P.base_hrs_sigma;          goto double_var;
			case 'm': sweep_ptr = &P.start_balance;           goto double_var;
			case 'r': sweep_ptr = &P.rich_threshold;          goto double_var;
			case 't': sweep_ptr = &P.savings_target;          goto double_var;
			case 'f': sweep_ptr = &P.savings_target_factor;   goto double_var;

			double_var:;
				*sweep_ptr = strtod( optarg, 0 );
				break;

			case 'k': P.max_kids = strtol( optarg, 0, 0 ); break;
			case 'o': P.overspend_ok = strtol( optarg, 0, 0 ); break;
			case 'z': sweep_step = 0.01; break;

			case 'q':
				if( !sweep_ptr )
				{
					fprintf( stderr, "illegal sweep\n" );
					exit( -1 );
				}
				sweep_from[ sweep_top ] = *sweep_ptr;
				sweep_to[ sweep_top ] = strtod( optarg, 0 );
				sweep_var[ sweep_top ] = sweep_ptr;
				++sweep_top;
				sweep_ptr = 0;
				break;

			case 'x': output_file = optarg; break;

			default:
				fprintf( stderr, "USAGE: %s [options]\n", argv[ 0 ]);
				for( c = 0; long_options[ c ].name; c++ )
				{
					fprintf( stderr, "    --%s%s\n", long_options[ c ].name, long_options[ c ].has_arg ? "=<val>" : "" );
				}
				exit( -1 );

		}
	}

	if( sweep_top )
	{
		sweeper( &P, sweep_top - 1 );
		exit( 0 );
	}

	setup_start( &P, &S );

	timeplot_init( &P );
	S.daytots = &daytot_timeplot;

	for( i = 0; i < TOT_DAYS; i++ )
	{
		hours_round( &P, &S );
		phonecall_round( &P, &S );
		collect_stats( &S );
		S.days ++;
	}
	timeplot_finish();

	printf( "Natural demand: %.2f hrs per person per day\n", S.demand_hrs / ( NUM_COUPLES * TOT_DAYS ));
	printf( "Traded hours: %.2f hrs per person per day\n", S.traded_hrs / ( NUM_COUPLES * TOT_DAYS ));
	printf( "Percent demand fulfillment: %.2f\n", 100 * S.traded_hrs / S.demand_hrs );
	return( 0 );
}

/*
 * total traded hours gives some approx equivalent to GDP and 
 * sort of represents the "success" of the babysitting scheme
 * (in as much as no one can go out without a babysitter and
 * more babysitting means more people get to go out)
 *
 * Various state counters consider the number of people who
 * wanted to go out but didn't because of lack of funds or because
 * of inability to find a babysitter.
 */