I need to find a permutation vector such that a 250x250 matrix reaches a maximal cost (sum of the elements in the superior triangle). Consider the cost function as a black box for my question. My neighborhood strategies are:
exchange 2 elements of the permutation vector
make an insertion of an element of the permutation vector elsewhere
I also have 2 pivotation rules:
first improvement
best improvement
My problem is: the loop is too long to execute. I should have a 20x slower code for this exercise, and I don't see were I could optimize my loop.
To do so, I use 3 loops. Here is the pseudo code (transform() can be either exchange, or insert, so a function using 2 indices (k, l) of the permutation vector currentSolution):
currentSolution = randomSolution(randomSeed)
bestCost = cost(currentSolution)
reapeat:
for each element k in currentSolution:
for each element l in currentSolution:
currentSolution = transform(currentSolution, k, l)
newCost = cost(currentSolution)
if newCost > bestCost:
depending on pivotation rule, update parameters (e.g. best_k, bestCost, ...)
depending on pivotation rule, stop the loops or revert transformation
else:
revert transformation of currentSolution
end
end
if at least one tuple (k, l) improves currentSolution:
currentSolution = transform(currentSolution, best_k, best_l) if not already done
else:
Found local optimum. Stop the search.
end
Now, here is the C code for the best improvement (insertion strategy):
while (i < EARLY_STOPPING) {
i++;
best_k = -1; best_l = -1; breakOuterLoop = 0;
for (k = 0; k <= PSize - 1; k++) {
for (l = 0; l <= PSize - 1; l++) {
if (l == k) continue;
insertElement(currentSolution, k, l);
newCost = computeCost(currentSolution);
if (newCost > bestCost) {
bestCost = newCost; best_k = k; best_l = l;
}
undoInsertion(currentSolution, k, l); // Put it back in place anyway
}
if (breakOuterLoop) break;
}
if (best_k >= 0) insertElement(currentSolution, best_k, best_l); // found a better k, l
else { // local optimum
break;
}
if (i >= EARLY_STOPPING) {
break;
}
}
NB: I know I could improve computeCost, but I am looking for another solution.
Here are the implementations of insertElement and undoInsertion, though I don't think there is a problem here:
void insertElement(long int* solution, int k, int l) {
if (k < l) {
long int temp = solution[k];
for (int i = k; i < l; i++) {
solution[i] = solution[i + 1];
}
solution[l] = temp;
} else if (k > l) {
long int temp = solution[k];
for (int i = k; i > l; i--) {
solution[i] = solution[i - 1];
}
solution[l] = temp;
}
// No operation needed if k == l, as an element cannot be inserted into its current position
}
void undoInsertion(long int* solution, int k, int l) {
if (k < l) {
long int temp = solution[l];
for (int i = l; i > k; i--) {
solution[i] = solution[i - 1];
}
solution[k] = temp;
} else if (k > l) {
long int temp = solution[l];
for (int i = l; i < k; i++) {
solution[i] = solution[i + 1];
}
solution[k] = temp;
}
// No operation needed if k == l, as an element cannot be "un-inserted" from its current position
}
How can I make my execution faster?
Edit: Here is the Makefile for the ones wondering how I compiled this project:
CC=gcc
CFLAGS=-O3 -Wall
OBJECTS=src/instance.o src/main.o src/optimization.o src/timer.o src/utilities.o
.PHONY: clean
all: lop
lop: $(OBJECTS)
$(CC) $(CFLAGS) $(OBJECTS) -o lop
clean:
rm -f src/*~ src/*.o lop
Here is also the computeCost function:
long long int computeCost (long int *s ) {
int h,k;
long long int sum;
/* Diagonal value are not considered */
for (sum = 0, h = 0; h < PSize; h++ )
for ( k = h + 1; k < PSize; k++ )
sum += CostMat[s[h]][s[k]];
return(sum);
}
Second edit: I have been asked the full code. Here is main.c:
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <string.h>
#include "instance.h"
#include "utilities.h"
#include "timer.h"
#include "optimization.h"
#define EARLY_STOPPING 20000
// Constants for first improvement pruning
#define IS_FIRST_THRESHOLD_MUTABLE 1
#define MIN_FIRST_THRESHOLD 50
#define FIRST_THRESHOLD_FACTOR_REDUCTION 2
long int FIRST_THRESHOLD = 10000; // this value must stay mutable
// Global variables for options
char *FileName;
char message[100];
int verbose = 0; // Verbose flag
int pivotingRule = 0; // 0 for first-improvement, 1 for best-improvement
int neighborhoodStrategy = 0; // 0 for transpose, 1 for exchange, 2 for insert
int initialSolutionMethod = 0; // 0 for random permutation, 1 for CW heuristic
void readOpts(int argc, char **argv) {
char opt;
FileName = NULL;
char *pivotingRuleStr, *neighborhoodStrategyStr, *initialSolutionMethodStr;
while ((opt = getopt(argc, argv, "i:p:n:s:v")) != -1) {
switch (opt) {
case 'i': // Instance file
FileName = strdup(optarg);
break;
case 'p': // Pivoting rule
pivotingRule = atoi(optarg);
pivotingRuleStr = (pivotingRule == 0) ? "first-improvement" : "best-improvement";
break;
case 'n': // Neighborhood strategy
neighborhoodStrategy = atoi(optarg);
switch(neighborhoodStrategy) {
case 0: neighborhoodStrategyStr = "transpose"; break;
case 1: neighborhoodStrategyStr = "exchange"; break;
case 2: neighborhoodStrategyStr = "insert"; break;
default: neighborhoodStrategyStr = "unknown"; break;
}
break;
case 's': // Initial solution method
initialSolutionMethod = atoi(optarg);
initialSolutionMethodStr = (initialSolutionMethod == 0) ? "random permutation" : "CW heuristic";
break;
case 'v': // Verbose mode
verbose = 1;
break;
default:
fprintf(stderr, "Option -%c not recognized.\n", opt);
}
}
if (verbose) {
if (FileName != NULL) {
verbose_message("Input file selected");
}
sprintf(message, "Pivoting rule selected is %s", pivotingRuleStr);
verbose_message(message);
sprintf(message, "Neighborhood strategy selected is %s", neighborhoodStrategyStr);
verbose_message(message);
sprintf(message, "Initial solution method selected is %s", initialSolutionMethodStr);
verbose_message(message);
}
if (!FileName) {
fprintf(stderr, "No instance file provided (use -i <instance_name>). Exiting.\n");
exit(1);
}
}
int main (int argc, char **argv)
{
long int i,j, k, l, step;
long int *currentSolution;
int cost, newCost, best_k, best_l, bestCost, oldCost;
int breakOuterLoop = 0; // Flag to signal breaking out of the outer loop
/* Do not buffer output */
setbuf(stdout,NULL);
setbuf(stderr,NULL);
if (argc < 2) {
printf("No instance file provided (use -i <instance_name>). Exiting.\n");
exit(1);
}
/* Read parameters */
readOpts(argc, argv);
/* Read instance file */
CostMat = readInstance(FileName);
printf("Data have been read from instance file. Size of instance = %ld.\n\n", PSize);
/* initialize random number generator, deterministically based on instance.
* To do this we simply set the seed to the sum of elements in the matrix, so it is constant per-instance,
but (most likely) varies between instances */
Seed = (long int) 0;
for (i=0; i < PSize; ++i)
for (j=0; j < PSize; ++j)
Seed += (long int) CostMat[i][j];
printf("Seed used to initialize RNG: %ld.\n\n", Seed);
/* starts time measurement */
start_timers();
/* A solution is just a vector of int with the same size as the instance */
currentSolution = (long int *)malloc(PSize * sizeof(long int));
if (initialSolutionMethod == 0) {
verbose_message(" Creating random solution...");
createRandomSolution(currentSolution);
verbose_message(" Done.");
} else if (initialSolutionMethod == 1) {
verbose_message(" Creating CW heuristic based solution...");
/* Initialize the attractiveness array and a flag array to keep track of inserted elements */
long int *attractiveness = (long int *)malloc(PSize * sizeof(long int));
int *inserted = (int *)calloc(PSize, sizeof(int)); // calloc initializes to 0
for (i = 0; i < PSize; i++) {
attractiveness[i] = 0;
for (k=0; k < PSize; k++) {
for (j = k; j < PSize; j++) {
attractiveness[i] += CostMat[i][j];
}
}
}
for (int step = 0; step < PSize; step++) {
// Find the most attractive element that hasn't been inserted yet
long int maxAttr = -1;
int maxIdx = -1;
for (i = 0; i < PSize; i++) {
if (!inserted[i] && (maxAttr < attractiveness[i])) {
maxAttr = attractiveness[i];
maxIdx = i;
}
}
inserted[maxIdx] = 1; // inserted
currentSolution[step] = maxIdx;
/* No need to update attractiveness since it's only used to determine the initial insertion order */
}
free(attractiveness);
free(inserted);
verbose_message(" CW heuristic based solution created.");
// Compute and print the cost of the CW heuristic based initial solution
newCost = computeCost(currentSolution);
sprintf(message, "Cost of the CW heuristic based initial solution = %d\n", newCost);
verbose_message(message);
} else {
printf("Unknown method for initial solution (should be 0 or 1, got %i).", initialSolutionMethod);
}
/* Print solution */
printf("Initial solution:\n");
for (j=0; j < PSize; j++)
printf(" %ld", currentSolution[j]);
printf("\n");
/* Compute cost of solution and print it */
cost = computeCost(currentSolution);
bestCost = cost;
sprintf(message, "Cost of the initial solution = %d\n", cost);
verbose_message(message);
swapAdjacentElements(currentSolution, 56);
bestCost = cost + computeSwapCostDifference(currentSolution, 56, 57);
verbose_message(" Starting main loop...");
oldCost = cost;
i = 0;
if (neighborhoodStrategy == 0) { // transpose
if (pivotingRule) { // transpose best improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1; // -1 means no k improves the model
long long int costDifference;
for (k = 0; k < PSize - 1; k++) { // Note the change here to ensure we don't go out of bounds
// Predict the new cost without actually swapping first
costDifference = computeSwapCostDifference(currentSolution, k, k + 1);
printf("%i\n", costDifference);
if (costDifference > 0) {
sprintf(message, "Found a neighbor augmenting the cost by %lli, new cost is %lli",
costDifference, bestCost + costDifference);
verbose_message(message);
bestCost = oldCost + costDifference; // Update the best known cost
best_k = k; // Remember the position that leads to the improvement
}
// No need to swap back since we didn't actually swap; we calculated the cost difference directly
}
if (best_k >= 0) {
// Now perform the actual swap for the best improvement found
swapAdjacentElements(currentSolution, best_k);
sprintf(message, "Performing swap at position %d for improvement.", best_k);
verbose_message(message);
} else {
verbose_message("Found a local optimum!");
break; // Local optimum reached
}
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
} else { // transpose first improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1;
for (k = 0; k <= PSize - 1; k++) {
swapAdjacentElements(currentSolution, k);
newCost = computeCost(currentSolution);
if (newCost > bestCost + FIRST_THRESHOLD) {
sprintf(message, "Found a neighbor augmenting the cost by %d, new cost is %d",
newCost - bestCost, newCost);
verbose_message(message);
bestCost = newCost;
best_k = k;
break; // Don't put back in place and leave the k-loop
} else {
swapAdjacentElements(currentSolution, k); // Put it back in place if no improvement
}
}
// If best k found, don't do anything as we didn't revert the change earlier
if (best_k < 0) {
if (IS_FIRST_THRESHOLD_MUTABLE && FIRST_THRESHOLD >= MIN_FIRST_THRESHOLD) {
FIRST_THRESHOLD /= FIRST_THRESHOLD_FACTOR_REDUCTION;
sprintf(message, "FIRST_THRESHOLD reduced to %li.", FIRST_THRESHOLD); verbose_message(message);
}
else {
verbose_message("Found a local optimum!");
break;
}
}
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
}
}
else if (neighborhoodStrategy == 1) { // exchange
if (pivotingRule) { // exchange best improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1; best_l = -1; breakOuterLoop = 0;
for (k = 0; k <= PSize - 1; k++) {
for (l = 0; l <= PSize - 1; l++) {
if (l == k) continue;
swapElements(currentSolution, k, l);
newCost = computeCost(currentSolution);
if (newCost > bestCost) {
sprintf(message, "Found a neighbor augmenting the cost by %d, new cost is %d",
newCost - bestCost, newCost);
verbose_message(message);
bestCost = newCost; best_k = k; best_l = l;
}
swapElements(currentSolution, k, l); // Put it back in place anyway
}
if (breakOuterLoop) break;
}
if (best_k >= 0) swapElements(currentSolution, best_k, best_l); // found a better k, l
else {
verbose_message("Found a local optimum!");
break;
} // local optimum
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
} else { // exchange first improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1; best_l = -1;
breakOuterLoop = 0;
for (k = 0; k <= PSize - 1; k++) {
best_l = -1;
for (l = 0; l <= PSize - 1; l++) {
if (l == k) continue; // Skip comparing an element with itself
swapElements(currentSolution, k, l);
newCost = computeCost(currentSolution);
if (newCost > bestCost + FIRST_THRESHOLD) {
sprintf(message, "Found a neighbor augmenting the cost by %d, new cost is %d",
newCost - bestCost, newCost);
verbose_message(message);
bestCost = newCost; best_k = k; best_l = l;
breakOuterLoop = 1; // Signal to break out of the outer loop
break; // Exit the inner l-loop
} else {
swapElements(currentSolution, k, l); // Put it back in place if no improvement
}
}
if (breakOuterLoop) break; // Check the flag and exit the outer k-loop if signaled
}
// If best k found, don't do anything as we didn't revert the change earlier
if (best_k < 0) {
if (IS_FIRST_THRESHOLD_MUTABLE && FIRST_THRESHOLD >= MIN_FIRST_THRESHOLD) {
FIRST_THRESHOLD /= FIRST_THRESHOLD_FACTOR_REDUCTION;
sprintf(message, "FIRST_THRESHOLD reduced to %li.", FIRST_THRESHOLD); verbose_message(message);
}
else {
verbose_message("Found a local optimum!");
break;
}
}
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
}
}
else if (neighborhoodStrategy == 2) { // insert
if (pivotingRule) { // insert best improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1; best_l = -1; breakOuterLoop = 0;
for (k = 0; k <= PSize - 1; k++) {
for (l = 0; l <= PSize - 1; l++) {
if (l == k) continue;
insertElement(currentSolution, k, l);
newCost = computeCost(currentSolution);
if (newCost > bestCost) {
sprintf(message, "Found a neighbor augmenting the cost by %d, new cost is %d",
newCost - bestCost, newCost);
verbose_message(message);
bestCost = newCost; best_k = k; best_l = l;
}
undoInsertion(currentSolution, k, l); // Put it back in place anyway
}
if (breakOuterLoop) break;
}
if (best_k >= 0) insertElement(currentSolution, best_k, best_l); // found a better k, l
else {
verbose_message("Found a local optimum!");
break;
} // local optimum
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
} else { // insert first improvement
while (i < EARLY_STOPPING) {
sprintf(message, "Starting iteration n°%li...", i);
verbose_message(message);
i++;
best_k = -1;
breakOuterLoop = 0;
for (k = 0; k <= PSize - 1; k++) {
best_l = -1;
for (l = 0; l <= PSize - 1; l++) {
if (l == k) continue; // Skip comparing an element with itself
insertElement(currentSolution, k, l);
newCost = computeCost(currentSolution);
if (newCost > bestCost + FIRST_THRESHOLD) {
sprintf(message, "Found a neighbor augmenting the cost by %d, new cost is %d",
newCost - bestCost, newCost);
verbose_message(message);
bestCost = newCost; best_k = k; best_l = l;
breakOuterLoop = 1; // Signal to break out of the outer loop
break; // Exit the inner l-loop
} else {
undoInsertion(currentSolution, k, l); // Put it back in place if no improvement
}
}
if (breakOuterLoop) break; // Check the flag and exit the outer k-loop if signaled
}
// If best k found, don't do anything as we didn't revert the change earlier
if (best_k < 0) {
if (IS_FIRST_THRESHOLD_MUTABLE && FIRST_THRESHOLD >= MIN_FIRST_THRESHOLD) {
FIRST_THRESHOLD /= FIRST_THRESHOLD_FACTOR_REDUCTION;
sprintf(message, "FIRST_THRESHOLD reduced to %li.", FIRST_THRESHOLD); verbose_message(message);
}
else {
verbose_message("Found a local optimum!");
break;
}
}
if (i >= EARLY_STOPPING) {
sprintf(message, "EARLY STOPPING (reached %li iterations).", i);
verbose_message(message);
break;
}
}
}
}
verbose_message(" Done.");
printf("Solution after optimization:\n");
for (j=0; j < PSize; j++)
printf(" %ld", currentSolution[j]);
printf("\n");
/* Recompute cost of solution after the exchange move */
/* There are some more efficient way to do this, instead of recomputing everything... */
newCost = computeCost(currentSolution);
sprintf(message, "Cost of the new solution = %i\n", newCost);
verbose_message(message);
cost = oldCost;
if (newCost == cost) {
verbose_message("Second solution is as good as first one");
}
else if (newCost > cost) {
sprintf(message, "Second solution is better than first one by %.2f%%.", ((float)(newCost - cost) / cost) * 100.0);
verbose_message(message); }
else {
verbose_message("Second solution is worse than first one");
}
printf("Time elapsed since we started the timer: %g\n\n", elapsed_time(VIRTUAL));
verbose_message("'main' completed successfully.");
printf("%s, %f, %d, %d, %.4f\n", FileName, elapsed_time(VIRTUAL), cost, newCost, ((float)(newCost - cost) / cost));
return 0;
}