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}
University of Michigan Program for the Study of Complex Systems
Contact http@maria.physics.lsa.umich.edu.
Revised November 4, 1996.