#!/usr/bin/perl -w
#
# popmodel:
# Very simple population model.
#
# Copyright (c) 2003 Chris Lightfoot. All rights reserved.
# Email: chris@ex-parrot.com; WWW: http://www.ex-parrot.com/~chris/
#

my $rcsid = ''; $rcsid .= '$Id:$';

use strict;

# initial age structure from
# http://www.national-statistics.gov.uk/census2001/pop2001/united_kingdom_ages.asp
my @agestructure = qw(
660080
680725
699913
712460
733291
721433
726866
747704
757129
785028
791067
774646
770701
783855
760340
763913
758448
726698
712268
702572
745538
739260
706745
673770
680838
700076
730325
760105
814122
862487
878633
869849
902531
908253
934319
930987
942077
933667
920581
898498
874182
840271
828071
818896
790160
762294
744666
753123
741433
734448
739798
765279
789112
852869
893379
689196
710966
689867
659506
589326
559590
591509
586244
579076
563529
547286
530226
507487
505819
506025
503581
486241
466129
443669
439611
426751
405387
389535
372571
372685
368167
337780
221508
188740
197352
184524
177304
153053
130373
107533);

my $y = 0;

sub sum ($$$) {
    my ($ary, $n1, $n2) = @_;
    my $S = 0;
    for (my $i = $n1; $i <= $n2; ++$i) {
        $S += $ary->[$i];
    }
    return $S;
}

sub update_pyramid ($) {
    my $py = shift;
    my @newpy;
    my $i;

    # 1. shift everyone to the right.
    $newpy[0] = 0;
    for ($i = 1; $i <= @$py; ++$i) {
        $newpy[$i] = $py->[$i - 1];
    }

    # 2. apply death rates; from http://www.statistics.gov.uk/STATBASE/ssdataset.asp?vlnk=6878
    my @deathrates = (
        0.25,       # 1--4
        0.1,        # 5--9
        0.15,
        0.4,
        0.55,
        0.6,
        0.8,
        1.05,
        1.55,
        2.55,
        3.7,
        6.55,
        10.25,
        16.7,
        28.7,
        49.25,
        84.6,
        174.4       # 85 and over
        );
    
    for ($i = 0; $i < @deathrates; ++$i) {
        my ($a, $b) = (1, 4);
        if ($i >= 1) {
            $a = 5 * $i;
            $b = $a + 4;
            $b = 200 if ($a == 85)
        }
        my $r = $deathrates[$i] / 1000.;
        for (my $j = $a; $j <= $b; ++$j) {
            $newpy[$j] ||= 0;
            $newpy[$j] *= 1 - $r;
            $newpy[$j] = 0 if ($newpy[$j] < 0);
        }
    }

    # 3. apply birth rate; from http://www.statistics.gov.uk/STATBASE/ssdataset.asp?vlnk=6893
    # nb / 2 since these are per woman.
    $newpy[0] = (27. / 1000.) * sum($py, 0, 19) / 2
                + (69.7 / 1000.) * sum($py, 20, 24) / 2.
                + (91.7 / 1000.) * sum($py, 25, 29) / 2.
                + (90.1 / 1000.) * sum($py, 30, 34) / 2.
                + (43.3 / 1000.) * sum($py, 35, 39) / 2.
                + (9.1 / 1000.) * sum($py, 40, 44) / 2.;
    # deaths per live births.
    $newpy[0] *= (1000 - 5.3) / 1000.;

    my $birthfactor = $ARGV[1];
    $birthfactor ||= 1;

    $newpy[0] *= $birthfactor;

    # 4. now add immigrants; from http://www.statistics.gov.uk/STATBASE/DatasetType.asp?vlnk=6914
    my @migration = (
        4000,       # 0 -- 14
        89000,      # 15 -- 24
        82000,      # 25 -- 44
        9000        # 45 and over
        );

    my $migrationfactor = $ARGV[0];
    $migrationfactor ||= 0;

    for ($i = 0; $i < 15; ++$i) {
        $newpy[$i] += $migrationfactor * $migration[0] / 14.;
    }
    for ($i = 15; $i < 25; ++$i) {
        $newpy[$i] += $migrationfactor * $migration[1] / 10.;
    }
    for ($i = 25; $i < 45; ++$i) {
        $newpy[$i] += $migrationfactor * $migration[2] / 20.;
    }
    # Assume no immigrants of age >65.
    for ($i = 45; $i < 65; ++$i) {
        $newpy[$i] += $migrationfactor * $migration[3] / 20.;
    }

    return \@newpy;
}

my $Y = 0;

my $a = \@agestructure;

my $worstyear = 0;
my $minratio = 5;
my $lastratio = 0;
my $firstratio = 0;

for ($Y = 0; $Y < 100; ++$Y) {
    if (exists($ENV{POPMODEL_PYRAMIDS})) {
        print join("\n", @$a);
        print "\n\n\n";
    }
    my $b = update_pyramid($a);
    $a = [@$b];

    my ($kids, $workers, $pensioners) = (sum($a, 0, 17), sum($a, 18, 64), sum($a, 65, 100));
    
    if (!exists($ENV{POPMODEL_PYRAMIDS})) {
        printf "%f %f %f %f %f\n", $kids + $workers + $pensioners, $kids, $workers, $pensioners, ($workers / $pensioners);
    }

    if ($Y > 10 && $Y < 40 && ($workers / $pensioners) < $minratio) {
        $minratio = $workers / $pensioners;
        $worstyear = $Y;
    }

    $lastratio = $workers / $pensioners;
    $firstratio = $workers / $pensioners if ($Y == 0);
}

print STDERR "$worstyear $minratio $lastratio $firstratio\n";