unit alliance_sim_calculate_unit;

interface

    uses
        alliance_sim_type_unit;

    procedure calculate_all_alliances (propensities: propensity_type; var potential_alliances: potential_alliance_type; top_alliance: longint; var global_index: longint; var tied_optima: integer; var permuted_propensities: propensity_type; var index_array: permuted_index_type; var tied_optima_array: tied_opt_list_type; var frustration_array: frustration_array_type; var optimum_array: optimum_array_type; var num_optima: integer; var size: size_type);

                            {------------------------------}

implementation
    procedure calculate_all_alliances (propensities: propensity_type; var potential_alliances: potential_alliance_type; top_alliance: longint; var global_index: longint; var tied_optima: integer; var permuted_propensities: propensity_type; var index_array: permuted_index_type; var tied_optima_array: tied_opt_list_type; var frustration_array: frustration_array_type; var optimum_array: optimum_array_type; var num_optima: integer; var size: size_type);

        type
            energy_pointer_type = ^blank_energy_rec;

            blank_energy_rec = record
                    value: real;
                    count: integer;
                    next: energy_pointer_type;
                    prev: energy_pointer_type;
                end;

        var
            country_i, country_j: integer;
            energy, current_energy, start_energy: real;
            current_bit, start_index: integer;

            global_energy: real;

            temp_current_alliance, current_alliance, alliance_config: longint;
            loop_count, x, y: longint;
            complement_alliance: longint;

                                                        {These are for permuting prop matrix}
            index: longint;
            new_matrix_i, new_matrix_j: longint;

                                                        {These are for frustration calculation}
            current_basin: longint;
            country, other_country: integer;
            temp_opt: opt_record;

            datetime: datetimerec;

                                        {-----------------------}

        function dist (alliance1, alliance2: boolean): integer;     {if in same alliance number, will get the "in" distance}
{                   if in diff, will get "between" distance.  The boolean is true if a bit in the rep is set, false if 0.}

                {If this distance is changed, change the message that prints out at the top of the output.}

        begin
            if (alliance1 = alliance2) then
                dist := 0
            else
                dist := 1;
        end;

                                        {-----------------------}

        procedure full_recursive_energy (var potential_alliances: potential_alliance_type; current_bit: integer; in_energy: real; current_index: longint);
            {if on last bit, saves this energy, otherwise sets the next bit to each of two values and calls self again.}
            {comes in with an index to a rep., and energy for that rep calculated up to but not including # current_bit. }
            {   Calculation currently does i <> j,  i.e. full matrix but no diagonal.  }

            var
                with_bit, next_bit: integer;
                current_energy: real;
                new_index: longint;

        begin
            {First, construct additional part of energy due to "current_bit" with all other bits.}
            {This calculation works directly with [i.e. loops on] the bit #s.  It could have been written using country #s}
            {   and the "bit" function I wrote, but it wasn't.  This works.  It saves one step of function calls this way.}

            current_energy := in_energy;

                {use the first calc line alone for the i < j matrix.  Current_bit makes i <= j,  + 1makes i < j }
            for with_bit := (current_bit + 1) to (num_countries - 1) do
                begin                       {for any point, this will do current with the _previous_ positions (that are fixed now)}
                                    {   that is, with the higher end bits.  }
                    current_energy := current_energy + propensities[country_from_bit(current_bit), country_from_bit(with_bit)] * dist(btst(current_index, current_bit), btst(current_index, with_bit));
                {to do asymmetric (full) matrix, the following line adds the other side of the diagonal from prop matrix}
                    current_energy := current_energy + propensities[country_from_bit(with_bit), country_from_bit(current_bit)] * dist(btst(current_index, current_bit), btst(current_index, with_bit));
                end;

                {To do the diagonal, necessary for true 1..n, 1..n, although it won't affect optima, add the following:  }
        {current_energy := current_energy + propensities[country_from_bit(current_bit), country_from_bit(current_bit)] * dist(btst(current_index, current_bit), btst(current_index, current_bit));}
                {If the diagonal is added or other energy changes are made, check frustration section to make sure it matches}

            if current_bit = 0 then
                begin               {have worked from high through low bit, and added low-bit energy contrib.}
                            {so make energy of this index be the  energy of that config}
                    potential_alliances^[current_index].energy := current_energy;

                            {can also set this energy into the complement alliance so don't have to recalc it}
                    complement_alliance := bitxor(current_index, top_alliance);
                    potential_alliances^[complement_alliance].energy := potential_alliances^[current_index].energy;
                end

            else                    {current_bit <> #0, so need to set next to 0, 1  and continue calc}
                begin
                    next_bit := current_bit - 1;
                    new_index := current_index;

                    BCLR(new_index, next_bit);              {set to 0; this command should be redundant.}
                    full_recursive_energy(potential_alliances, next_bit, current_energy, new_index);

                    BSET(new_index, next_bit);              {set to 1}
                    full_recursive_energy(potential_alliances, next_bit, current_energy, new_index);

                end;

        end;                {proc full_recursive_energy}

            { ---------------------------------------  }

        procedure local_find (var potential_alliances: potential_alliance_type; index_in: longint; levels_deep: integer);
                {recursive proc to find the local optimum of an alliance in index_in}
            var
                best_adj: longint;
        begin
            if potential_alliances^[index_in].local_opt <> max_alliances_plus1 then
                begin
                        {DONE - HAD FOUND A LOCAL IN A PREVIOUS RECURSIVE SEARCH}
    {Dont assign anything new to potential.local_opt}
                { old statment when was a func:  local_find := potential_alliances^[index_in].local_opt;}
                end
            else                {have to look for the best local...}
                begin
                    best_adj := best_neighbor(index_in);
                    if best_adj = index_in then
                        begin                               {if best is yourself, then at a local opt.}
                            potential_alliances^[index_in].local_opt := best_adj;
                        end
                    else                                    {best must be searched from}
                        begin
                        {Call the procedure again to find the optimum of point best_adj; store that in potential array}
                            levels_deep := levels_deep + 1;
                            if levels_deep > num_countries then
                                begin
                                    writeln('Error in recursive local_find procedure -- calling deeper than num_countries');
                                end;
                            local_find(potential_alliances, best_adj, levels_deep);
                        {When the optimum of the best adjacent point comes out, that is also the local max of the current point}
                            potential_alliances^[index_in].local_opt := potential_alliances^[best_adj].local_opt;
                        end;
                end;
        end;                        {function local_find}

            { ---------------------------------------  }

        procedure add_one_count_to_optimum_array (an_index: longint; var optimum_array: optimum_array_type; num_optima: integer);
            {procedure takes one occurence of an_index as a local max and adds it on the count for that optimum}
            var
                location: integer;
        begin
                    {first figure out which optimum in the optimum_array  an_index really is.}
                    {  Note - this search could be sped up, but given that there are usually few optima, it's probably }
                    {     not a big problem, as it will usually do only a couple of steps.  }
            location := 0;
            repeat
                location := location + 1;
            until (an_index = optimum_array[location].index) or (location >= num_optima) or (location >= max_optima);
            if ((location >= max_optima) or (location >= num_optima)) and (an_index <> optimum_array[location].index) then
                begin
    {writeln('Error in proc  add_one_count_to_optimum_array -- location exceeded number of optima in list');}
    {This used to write above message, but since this (loc >= num) will happen anytime more than num_optima are }
    {  found, it no longer writes anything.  The message that more than the acceptable number of optima were found is }
    {   printed from the main calculate procedure.  }
                end
            else                        {OK - have position in list}
                begin
                    optimum_array[location].basin_size := optimum_array[location].basin_size + 1;
                end;
        end;                        {proc add one count}

            { ---------------------------------------  }

    begin               {main proc calc all alliances}
        gettime(datetime);
        writeln('Starting energy calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);

{First calc all alliance energies}
        for alliance_config := 0 to (top_alliance) do
            begin
                potential_alliances^[alliance_config].energy := 0;
                potential_alliances^[alliance_config].local_opt := max_alliances_plus1;
            end;

        {Previous energy calc replaced by this:  }
        current_bit := num_countries - 1;           {Start with high end bit for the number of countries we have}
        start_energy := 0;                                  {Start with a 0 energy set to add on to}
        start_index := 0;                                   {Starting with alliance index 0}

        full_recursive_energy(potential_alliances, current_bit, start_energy, start_index);
    {calling like this will figure energy of all non-complement alliances, then put into complements}


{Now calculate best adjacent point and follow gradients to get local optima}
        gettime(datetime);
        writeln('Starting local optimum calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);
        for alliance_config := 0 to top_alliance do
            begin
                local_find(potential_alliances, alliance_config, 0);
            end;
                    {Now have calculated local maxima by following gradients.}
        {Local find exits with the array of all ties and ranks (broken ties) still intact in first_tie_value}

                    {Now put the optima on the main optimum list.  }
        gettime(datetime);
        writeln('Starting optimum count operation at ', datetime.hour : 2, ': ', datetime.minute : 2);
        num_optima := 0;
                    {if a point is an optima, then  put it on the list}
        for current_alliance := 0 to top_alliance do
            if potential_alliances^[current_alliance].local_opt = current_alliance then
                if num_optima < max_optima then
                    begin                       {its an optimum, so add the point onto the optimum_array}
                        num_optima := num_optima + 1;
                        optimum_array[num_optima].index := current_alliance;
                    end
                else                    {num_optima >= max_optima, so adding one more would go over list.}
                    begin
                        if num_optima = max_optima then             {write message once}
                            begin
                                writeln('More than ', max_optima : 3, ' basins were seen.  The program can only handle ', max_optima : 3, '.  Others will be ignored.');
                                writeln;
                            end;
                        num_optima := num_optima + 1;
                    end;
        if odd(num_optima) then
            writeln('ERROR in program - counted an odd number of optima -- should be EVEN');

                {At end of above, num_optima is count in both base and complement set of how many configs}
                {  have basin size > 0.  So, it is number of optima in both halves.  This should be an even number.}
                {Optimum Array  will have num_optima entries in it, in positions 1 to num_optima.}
                {    Note added 11/6/90}
                {  Note remains true after modifications of 9/91. }


                {Now sort   __the list of basins__   by energy.  }
                {Previously sorted each half separately, since assumed complements.  But they don't have to be.}
                {  So, just sort whole list by rank, and just don't output the complements later if there are any there.}

        for x := 1 to min(num_optima, max_optima) - 1 do                {all}
            for y := x + 1 to min(num_optima, max_optima) do
                begin
                    if potential_alliances^[optimum_array[y].index].energy < potential_alliances^[optimum_array[x].index].energy then
                        begin                   {switch}
                            temp_opt := optimum_array[x];
                            optimum_array[x] := optimum_array[y];
                            optimum_array[y] := temp_opt;
                        end;
                end;


                    {Now calculate size of basins to go with the list constructed above.}
        gettime(datetime);
        writeln('Starting basin size calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);
        for current_alliance := 0 to top_alliance do
            add_one_count_to_optimum_array(potential_alliances^[current_alliance].local_opt, optimum_array, num_optima);



                                {Now, figure out global optimum and ties for it.}
        gettime(datetime);
        writeln('Starting global optimum calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);

        global_energy := potential_alliances^[0].energy;
        global_index := 0;
        for current_alliance := 1 to top_alliance do
            if potential_alliances^[current_alliance].energy < global_energy then
                {Since using >, this will keep the non-complement record as the global optimum.}
                begin
                    global_energy := potential_alliances^[current_alliance].energy;
                    global_index := current_alliance;
                end;


          {Now count any ties for the high (global optimum) value.  If there are any, they might be on the tie_value list.}
{           However, they also might not be, b/c values are only put on there if a point tries to go to both of them.}
{           So, to make sure I have all the ties, go through the full alliance list one more time.  }
        gettime(datetime);
        writeln('Starting global ties calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);
        tied_optima := 0;               {Initialize to 0, since going through whole list.}
                                            {should always get at least one more increment from being tied with the complement}
        for current_alliance := 0 to top_alliance do
            begin
                if potential_alliances^[current_alliance].energy = global_energy then
                    if tied_optima < max_tied_optima then
                        begin
                            tied_optima := tied_optima + 1;
                            tied_optimum_array[tied_optima] := current_alliance;
                        end
                    else                        {Seen more than max;  keep counting, but don't save.}
                        begin
                            if tied_optima = max_tied_optima then           {write message once}
                                writeln('More ties for global optimum than allowed.  Only the first ', max_tied_optima : 3, ' have been stored.');
                            tied_optima := tied_optima + 1;
                        end;
            end;                {for current_alliance 0 to top}
        if odd(tied_optima) then
            writeln('ERROR in program - counted an odd number of optima -- should be EVEN');

                {From the above section, the number of tied_optimas is the count for both sets - basic and}
{                  Complements, and should be correct even when reaching max limit, or reaching top alliance.  }
{                  The number should go up in Increments of 2, since counting basic plus complements.  }


        gettime(datetime);
        writeln('Starting propensity permutation and frustration calculation at ', datetime.hour : 2, ': ', datetime.minute : 2);

                {Now calculate permuted propensity matrix}
                {first build an index of the proper order of nations in the new propensity matrix}
        index := 1;
        for x := 1 to num_countries do
            if btst(global_index, bit(x)) = false then
                begin
                    index_array[index] := x;
                    index := index + 1;
                end;
        for x := 1 to num_countries do
            if btst(global_index, bit(x)) = true then
                begin
                    index_array[index] := x;
                    index := index + 1;
                end;
                        {now use index array to rearrange elements of the original prop matrix}
        for new_matrix_i := 1 to num_countries do
            for new_matrix_j := 1 to num_countries do
                begin
                    permuted_propensities[new_matrix_i, new_matrix_j] := propensities[index_array[new_matrix_i], index_array[new_matrix_j]];
                end;



                    {Now calculate frustrations of all countries for all optima}

                            {Compute actual "payoff" at each optimum to each country, and put in frustration matrix.}
{Frustration for country i is defined as sum over all j of size j * raw_prop(i,j) * distance(i,j)}
        for country := 1 to num_countries do
            begin
                for current_basin := 1 to min(num_optima, max_optima) do
                    begin
                        frustration_array[country, current_basin] := 0;
                        for other_country := 1 to num_countries do                  {calc contrib of country "other" to actual payoff}
                            if other_country <> country then                {don't do calc of the country w/itself, since ignore diagonal before}
                                        {add raw_prop * distance(me,other) *size (other)   for the config of this optimum}
                                frustration_array[country, current_basin] := frustration_array[country, current_basin] + ((propensities[country, other_country] / size[country]) / size[other_country]) * size[other_country] * dist(btst(optimum_array[current_basin].index, bit(country)), btst(optimum_array[current_basin].index, bit(other_country)));
                    end;                            {for current basin}
            end;


                                    {Additionally, make frustration [0] be frust with input alliance}
        if have_starting_alliance then
            for country := 1 to num_countries do
                begin
                    frustration_array[country, 0] := 0;
                    for other_country := 1 to num_countries do                      {calc contrib of country "other" to actual payoff}
                        if other_country <> country then                                {don't do calc of the country itself}
                            frustration_array[country, 0] := frustration_array[country, 0] + ((propensities[country, other_country] / size[country]) / size[other_country]) * size[other_country] * dist(btst(starting_alliance.index, bit(country)), btst(starting_alliance.index, bit(other_country)));
                                            {add prop * distance in the optimum}
                end;                            {for country}



    end;                    {procedure calc all alliances}


end.                {unit calculate}