Java tutorial
/* * SymmetricReplicatorDynamics.java * * Copyright (C) 2006-2009 Patrick R. Jordan * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package egat.replicatordynamics; import egat.game.*; import egat.minform.SymmetricRationalizableFinder; import egat.minform.LpSolveSymmetricRationalizableFinder; import java.util.Set; import java.util.Arrays; import java.io.PrintStream; import org.apache.commons.math.random.RandomDataImpl; import org.apache.commons.math.random.RandomData; /** * @author Patrick Jordan */ public class SymmetricConstrainedFeasibleAmoebaSearch { private double tolerance; private int maxIteration; private PrintStream printStream; private RandomData randomData; public SymmetricConstrainedFeasibleAmoebaSearch(double tolerance, int maxIteration, PrintStream printStream, SymmetricRationalizableFinder rationalizableFinder) { this.tolerance = tolerance; this.maxIteration = maxIteration; this.printStream = printStream; this.randomData = new RandomDataImpl(); } public SymmetricConstrainedFeasibleAmoebaSearch(double tolerance, int maxIteration, PrintStream printStream) { this(tolerance, maxIteration, printStream, new LpSolveSymmetricRationalizableFinder()); } public SymmetricConstrainedFeasibleAmoebaSearch(double tolerance, int maxIteration) { this(tolerance, maxIteration, null); } public Strategy run(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions) { return run(game, initialStrategy, restrictedActions, false); } public StrategySearchResult runForResult(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions) { return runForResult(game, initialStrategy, restrictedActions, false); } public Strategy run(Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize, Action[] actions, PayoffMatrix pm) { boolean[] mask = new boolean[actions.length]; int restrictedCount = 0; for (int i = 0; i < actions.length; i++) { mask[i] = restrictedActions.contains(actions[i]); if (mask[i]) { restrictedCount++; } } int[] restricted = new int[restrictedCount]; for (int i = 0, restrictedIndex = 0; i < actions.length; i++) { if (mask[i]) { restricted[restrictedIndex++] = i; } } // Normalize to a uniform distribution double[] distribution = new double[actions.length]; double[] bestDist = new double[actions.length]; if (initialStrategy == null) { double sum = 0.0; for (int i = 0; i < bestDist.length; i++) { if (mask[i]) { bestDist[i] = Math.random(); sum += bestDist[i]; } } for (int i = 0; i < bestDist.length; i++) { if (mask[i]) { bestDist[i] /= sum; } } } else { for (int i = 0; i < actions.length; i++) { if (mask[i]) { bestDist[i] = initialStrategy.getProbability(actions[i]).doubleValue(); } } } double bestRegret = pm.regret(bestDist); // Initialize current strategy Strategy currentStrategy = buildStrategy(actions, bestDist, restricted); fireUpdatedStrategy(currentStrategy, 0, Double.NaN); int SIMPLICES = restricted.length; double[][] X = new double[SIMPLICES][SIMPLICES - 1]; double[] centroid = new double[SIMPLICES - 1]; double[] reflect = new double[SIMPLICES - 1]; double[] expand = new double[SIMPLICES - 1]; double[] contract = new double[SIMPLICES - 1]; double[] V = new double[SIMPLICES]; double[] cur = distribution.clone(); if (!randomize) { for (int r = 0, n = SIMPLICES - 1; r < n; r++) { X[0][r] = 0.0; } V[0] = calculateObjective(pm, X[0], restricted, distribution); for (int s = 1; s < SIMPLICES; s++) { System.arraycopy(X[0], 0, X[s], 0, SIMPLICES - 1); X[s][s - 1] = 1.0; V[s] = calculateObjective(pm, X[s], restricted, distribution); } } else { for (int s = 0; s < SIMPLICES; s++) { generateSimplex(X[s]); V[s] = calculateObjective(pm, X[s], restricted, distribution); } } int highest = -1; int lowest = -1; int second = -1; for (int s = 0; s < SIMPLICES; s++) { if (lowest < 0 || V[s] < V[lowest]) { lowest = s; } if (second < 0 || V[s] > V[second]) { if (highest < 0 || V[s] > V[highest]) { second = highest; highest = s; } else { second = s; } } } double curRegret = Double.POSITIVE_INFINITY; loop: for (int iteration = 1; SIMPLICES > 1; iteration++) { generateCentroid(centroid, X, SIMPLICES, highest); generateFeasibleTowards(reflect, centroid, X[highest], 1.0, pm, restricted, distribution); double reflectV = calculateObjective(pm, reflect, restricted, distribution); // Reflect if (V[lowest] <= reflectV && reflectV < V[second]) { replaceHighest(reflect, reflectV, X, V, highest); // Expand } else if (reflectV < V[lowest]) { generateFeasibleTowards(expand, centroid, X[highest], 2.0, pm, restricted, distribution); double expandV = calculateObjective(pm, expand, restricted, distribution); if (expandV < reflectV) { replaceHighest(expand, expandV, X, V, highest); } else { replaceHighest(reflect, reflectV, X, V, highest); } // Contract } else { generateAway(contract, X[highest], centroid, 0.5); double contractV = calculateObjective(pm, contract, restricted, distribution); if (contractV < V[highest]) { replaceHighest(contract, contractV, X, V, highest); // Reduction } else { for (int s = 0; s < SIMPLICES; s++) { if (s != lowest) { generateAway(X[s], X[lowest], X[s], 0.5); V[s] = calculateObjective(pm, X[s], restricted, distribution); } } } } highest = -1; lowest = -1; second = -1; for (int s = 0; s < SIMPLICES; s++) { if (lowest < 0 || V[s] < V[lowest]) { lowest = s; } if (second < 0 || V[s] > V[second]) { if (highest < 0 || V[s] > V[highest]) { second = highest; highest = s; } else { second = s; } } } Arrays.fill(cur, 0.0); double sum = 0.0; for (int r = 0; r < restricted.length - 1; r++) { cur[restricted[r]] = Math.max(0, X[lowest][r]); sum += cur[restricted[r]]; } cur[restricted[restricted.length - 1]] = Math.max(0.0, 1 - sum); sum += cur[restricted[restricted.length - 1]]; for (int r = 0; r < restricted.length; r++) { cur[restricted[r]] /= sum; } curRegret = pm.regret(cur); if (curRegret < bestRegret) { System.arraycopy(cur, 0, bestDist, 0, cur.length); bestRegret = curRegret; } // Build the new strategy currentStrategy = buildStrategy(actions, bestDist, restricted); fireUpdatedStrategy(currentStrategy, iteration, bestRegret); // Calculate the norm double norm = V[highest] - V[lowest]; // Check termination condition if (terminate(norm, iteration)) break; } return currentStrategy; } public StrategySearchResult runForResult(Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize, Action[] actions, PayoffMatrix pm) { boolean[] mask = new boolean[actions.length]; int restrictedCount = 0; for (int i = 0; i < actions.length; i++) { mask[i] = restrictedActions.contains(actions[i]); if (mask[i]) { restrictedCount++; } } int[] restricted = new int[restrictedCount]; for (int i = 0, restrictedIndex = 0; i < actions.length; i++) { if (mask[i]) { restricted[restrictedIndex++] = i; } } // Normalize to a uniform distribution double[] distribution = new double[actions.length]; double[] bestDist = new double[actions.length]; if (initialStrategy == null) { double sum = 0.0; for (int i = 0; i < bestDist.length; i++) { if (mask[i]) { bestDist[i] = Math.random(); sum += bestDist[i]; } } for (int i = 0; i < bestDist.length; i++) { if (mask[i]) { bestDist[i] /= sum; } } } else { for (int i = 0; i < actions.length; i++) { if (mask[i]) { bestDist[i] = initialStrategy.getProbability(actions[i]).doubleValue(); } } } double bestRegret = pm.regret(bestDist); // Initialize current strategy Strategy currentStrategy = buildStrategy(actions, bestDist, restricted); fireUpdatedStrategy(currentStrategy, 0, Double.NaN); int SIMPLICES = restricted.length; double[][] X = new double[SIMPLICES][SIMPLICES - 1]; double[] centroid = new double[SIMPLICES - 1]; double[] reflect = new double[SIMPLICES - 1]; double[] expand = new double[SIMPLICES - 1]; double[] contract = new double[SIMPLICES - 1]; double[] V = new double[SIMPLICES]; double[] cur = distribution.clone(); if (!randomize) { for (int r = 0, n = SIMPLICES - 1; r < n; r++) { X[0][r] = 0.0; } V[0] = calculateObjective(pm, X[0], restricted, distribution); for (int s = 1; s < SIMPLICES; s++) { System.arraycopy(X[0], 0, X[s], 0, SIMPLICES - 1); X[s][s - 1] = 1.0; V[s] = calculateObjective(pm, X[s], restricted, distribution); } } else { for (int s = 0; s < SIMPLICES; s++) { generateSimplex(X[s]); V[s] = calculateObjective(pm, X[s], restricted, distribution); } } int highest = -1; int lowest = -1; int second = -1; for (int s = 0; s < SIMPLICES; s++) { if (lowest < 0 || V[s] < V[lowest]) { lowest = s; } if (second < 0 || V[s] > V[second]) { if (highest < 0 || V[s] > V[highest]) { second = highest; highest = s; } else { second = s; } } } double curRegret = Double.POSITIVE_INFINITY; int iteration = 0; double norm = Double.POSITIVE_INFINITY; loop: for (iteration = 1; SIMPLICES > 1; iteration++) { generateCentroid(centroid, X, SIMPLICES, highest); generateFeasibleTowards(reflect, centroid, X[highest], 1.0, pm, restricted, distribution); double reflectV = calculateObjective(pm, reflect, restricted, distribution); // Reflect if (V[lowest] <= reflectV && reflectV < V[second]) { replaceHighest(reflect, reflectV, X, V, highest); // Expand } else if (reflectV < V[lowest]) { generateFeasibleTowards(expand, centroid, X[highest], 2.0, pm, restricted, distribution); double expandV = calculateObjective(pm, expand, restricted, distribution); if (expandV < reflectV) { replaceHighest(expand, expandV, X, V, highest); } else { replaceHighest(reflect, reflectV, X, V, highest); } // Contract } else { generateAway(contract, X[highest], centroid, 0.5); double contractV = calculateObjective(pm, contract, restricted, distribution); if (contractV < V[highest]) { replaceHighest(contract, contractV, X, V, highest); // Reduction } else { for (int s = 0; s < SIMPLICES; s++) { if (s != lowest) { generateAway(X[s], X[lowest], X[s], 0.5); V[s] = calculateObjective(pm, X[s], restricted, distribution); } } } } highest = -1; lowest = -1; second = -1; for (int s = 0; s < SIMPLICES; s++) { if (lowest < 0 || V[s] < V[lowest]) { lowest = s; } if (second < 0 || V[s] > V[second]) { if (highest < 0 || V[s] > V[highest]) { second = highest; highest = s; } else { second = s; } } } Arrays.fill(cur, 0.0); double sum = 0.0; for (int r = 0; r < restricted.length - 1; r++) { cur[restricted[r]] = Math.max(0, X[lowest][r]); sum += cur[restricted[r]]; } cur[restricted[restricted.length - 1]] = Math.max(0.0, 1 - sum); sum += cur[restricted[restricted.length - 1]]; for (int r = 0; r < restricted.length; r++) { cur[restricted[r]] /= sum; } curRegret = pm.regret(cur); if (curRegret < bestRegret) { System.arraycopy(cur, 0, bestDist, 0, cur.length); bestRegret = curRegret; } // Build the new strategy currentStrategy = buildStrategy(actions, bestDist, restricted); fireUpdatedStrategy(currentStrategy, iteration, bestRegret); // Calculate the norm norm = V[highest] - V[lowest]; // Check termination condition if (terminate(norm, iteration)) break; } return new StrategySearchResult(currentStrategy, bestRegret, iteration, norm); } public Strategy run(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize, Player[] players, Action[] actions) { return run(initialStrategy, restrictedActions, randomize, actions, createPayoffMatrix(game, players, actions)); } public StrategySearchResult runForResult(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize, Player[] players, Action[] actions) { return runForResult(initialStrategy, restrictedActions, randomize, actions, createPayoffMatrix(game, players, actions)); } public Strategy run(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize) { Action[] actions = game.getActions().toArray(new Action[0]); Player[] players = game.players().toArray(new Player[0]); return run(initialStrategy, restrictedActions, randomize, actions, createPayoffMatrix(game, players, actions)); } public StrategySearchResult runForResult(SymmetricGame game, Strategy initialStrategy, Set<Action> restrictedActions, boolean randomize) { Action[] actions = game.getActions().toArray(new Action[0]); Player[] players = game.players().toArray(new Player[0]); return runForResult(initialStrategy, restrictedActions, randomize, actions, createPayoffMatrix(game, players, actions)); } private void replaceHighest(double[] src, double srcV, double X[][], double[] V, int highest) { System.arraycopy(src, 0, X[highest], 0, src.length); V[highest] = srcV; } private void generateCentroid(double[] centroid, double[][] X, int simplices, int highest) { for (int r = 0; r < simplices - 1; r++) { centroid[r] = 0.0; for (int s = 0; s < simplices; s++) { if (s != highest) { centroid[r] += X[s][r]; } } centroid[r] /= simplices - 1.0; } } private void generateFeasibleTowards(double[] newVertex, double[] first, double[] second, double scale, PayoffMatrix pm, int[] restricted, double[] distribution) { generateTowards(newVertex, first, second, scale); double V = calculateObjective(pm, newVertex, restricted, distribution); if (Double.isInfinite(V)) { double a1 = scale; double a2 = 0.0; double a; double aV; do { a = (a1 + a2) * 0.5; generateTowards(newVertex, first, second, a); aV = calculateObjective(pm, newVertex, restricted, distribution); if (Double.isInfinite(aV)) { a1 = a; } else { a2 = a; } } while (a1 - a2 > 1e-5); if (Double.isInfinite(aV)) { generateTowards(newVertex, first, second, a2); } } } private void generateTowards(double[] newVertex, double[] first, double[] second, double scale) { for (int r = 0; r < newVertex.length; r++) { newVertex[r] = first[r] + scale * (first[r] - second[r]); } } private void generateAway(double[] newVertex, double[] first, double[] second, double scale) { generateTowards(newVertex, first, second, -scale); } private double calculateObjective(PayoffMatrix pm, double[] cur, int[] restricted, double[] distribution) { double sum = 0.0; for (int r = 0, n = restricted.length - 1; r < n; r++) { if (cur[r] < 0) { return Double.POSITIVE_INFINITY; } sum += cur[r]; } if (sum > 1.0) { return Double.POSITIVE_INFINITY; } return calculateRegret(pm, cur, restricted, distribution); } private double calculateRegret(PayoffMatrix pm, double[] cur, int[] restricted, double[] distribution) { Arrays.fill(distribution, 0.0); double sum = 0.0; for (int r = 0; r < restricted.length - 1; r++) { distribution[restricted[r]] = cur[r]; sum += cur[r]; } distribution[restricted[restricted.length - 1]] = 1.0 - sum; return pm.regret(distribution); } private void fireUpdatedStrategy(Strategy currentStrategy, int iteration, double epsilon) { if (printStream != null) { printStream.println( String.format("(s: %s iteration: %d epsilon: %s", currentStrategy, iteration, epsilon)); } } private Strategy buildStrategy(Action[] actions, double[] distribution, int[] restricted) { Action[] maskedActions = new Action[restricted.length]; Number[] number = new Number[restricted.length]; for (int i = 0; i < restricted.length; i++) { maskedActions[i] = actions[restricted[i]]; number[i] = distribution[restricted[i]]; } return Games.createStrategy(maskedActions, number); } private boolean terminate(double norm, int iteration) { return (norm < tolerance || iteration > maxIteration); } private PayoffMatrix createPayoffMatrix(SymmetricGame game, Player[] players, Action[] actions) { Action[] playerActions = new Action[players.length]; PayoffMatrix pm = new PayoffMatrix(players.length, actions.length); int[] indices = new int[players.length]; for (int i = 0; i < pm.getTotalSize(); i++) { pm.expand(i, indices); for (int j = 0; j < players.length; j++) { playerActions[j] = actions[indices[j]]; } Outcome outcome = Games.createOutcome(players, playerActions); try { Payoff p = game.payoff(outcome); pm.setPayoff(i, p.getPayoff(players[0]).getValue()); } catch (NonexistentPayoffException e) { pm.setPayoff(i, -Double.MAX_VALUE); } } return pm; } private void generateSimplex(double[] X) { double sum = 0.0; for (int i = 0; i < X.length; i++) { X[i] = randomData.nextExponential(1.0); sum += X[i]; } sum += randomData.nextExponential(1.0); for (int i = 0; i < X.length; i++) { X[i] /= sum; } } }