2
\$\begingroup\$

I wrote this program to calculate a value from an energy function from a sequence of 1s and -1s. The problem here is to find the lowest possible fitness for a function (note that this is original function; in my implementation I pre calculate power of c c*c and put it in the thread so I avoid calling a method):

double EnergyFunction::evaluate(vector<short> field) {
  double e = 0.0;
  int l = field.size()
  for (int k = 1; k < l; k++){
    double c = 0.0;
    for (int i = 0, j = k; j < l; i++, j++) {
        c += field[i] * field[j];
    }
    e += c*c;
  }
  return l*l / ( e+e );
}

It uses multithreading and one thread printOut() for outputing data to standard output stream. In thread printOut() there is a mutex for outputting without a data race. With this I hoped to achieve no locking of other threads, where computation is done.

Concerns:

  1. Is this the right way to do it without locking other threads?

    The sequence is randomly mutated every loop iteration, so every iteration I am copying the sequence from the best sequence. That way I start every time with the current best sequence and mutate it to get a new fitness.

  2. Is copying the whole sequence every iteration necessary to get best solution every time, or is there any other way to get the same effect? I cannot think of it. Also bear in mind the algorithm isn't so important than performance. The point is still to achieve the lowest possible fitness but I want to make it as simple as possible and as fast as possible.

#include <iostream>
#include <bitset>
#include <atomic>
#include <random>
#include <chrono>
#include <array>
#include <thread>
#include <mutex>
#include <sstream>
#include <math.h>

static const bool lookupTable[2] = {false, true};
static const short lookupTableNum[2] = {-1, 1};
unsigned int lookupTableCPower[305] = {};
static const unsigned long iterations = 100000000;
static const short l = 305;
static const short lh = 153;

std::atomic<unsigned long> overallGeneration (1);
std::atomic<unsigned short> bestFitness (std::numeric_limits<unsigned short>::max());
std::array<std::atomic<bool>, lh> bestSequence;

std::random_device rd;
std::mt19937 mt(rd());

std::uniform_int_distribution<unsigned short> distDim(0.0, lh);
std::uniform_int_distribution<unsigned short> dist(0.0, 1.0);

std::mutex mtx;

void threadOptimize();
void printOut(std::chrono::high_resolution_clock::time_point begin);

int main(int argc, const char * argv[]) {
    for(int i = 0; i < 305; i++)
        lookupTableCPower[i] = i * i;

    int numberOfThreads = (int)std::thread::hardware_concurrency();
    std::thread *threads = new std::thread[numberOfThreads];

    for (int i = 0; i < numberOfThreads; i++)
        threads[i] = std::thread(threadOptimize);

    for (int i = 0; i < numberOfThreads; i++)
        threads[i].join();

    std::cout << "Program ended." << std::endl;

    delete [] threads;
    return 0;
}

void threadOptimize() {
    bool seq[l];
    for(short i = 0; i < lh; i++)
        seq[i] = lookupTable[dist(mt)];

    std::copy(seq, seq + lh, bestSequence.begin());
    auto begin = std::chrono::high_resolution_clock::now();

    while(overallGeneration.load(std::memory_order_acquire) < iterations) {
        for(int i = 0; i  < distDim(mt);++i)
            seq[distDim(mt)] ^= true;

        bool even = false;
        for(int i = 152, j = 154; i >= 0 && j < 305; i--,j++) {
            if(even)
                seq[j] = seq[i] ^ true;
            else
                seq[j] = seq[i];
            even ^= true;
        }

        unsigned int newFitness = 0;
        for (int k = 1; k < l; ++k) {
            short c = 0;
            for(int i = 0, j = k; j < l; ++i,++j) {
                c += lookupTableNum[(seq[i] ^ seq[j])];
            }
            newFitness += lookupTableCPower[(c >= 0) ? c : -c];
        }

        ++overallGeneration;
        if (newFitness < bestFitness.load(std::memory_order_acquire)) {
            bestFitness.store(newFitness,std::memory_order_release);
            std::copy(seq, seq+lh, bestSequence.begin());

            std::thread t1(printOut, begin);
            t1.detach();
        } else
            std::copy(bestSequence.begin(), bestSequence.end(), seq);
    }
}

void printOut(std::chrono::high_resolution_clock::time_point begin) {
    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end-begin).count();

    duration /= 1000;
    duration++;

    std::lock_guard<std::mutex> lock(mtx);
    std::cout << (overallGeneration.load(std::memory_order_acquire) / (duration)) << " " << bestFitness.load(std::memory_order_acquire) << " ";
    for(short i = 0; i < l; i++)
        std::cout << bestSequence[i].load(std::memory_order_acquire);
    std::cout << std::endl;
}

I am also trying to remove branch prediction and pre-calculate numbers for the energy function.

\$\endgroup\$
  • 1
    \$\begingroup\$ Feel free to post a follow-up with the revised code. Or, preferably, wait some more for other answers and thén post a follow-up with the edited code. \$\endgroup\$ – Mast Dec 3 '15 at 15:00
2
\$\begingroup\$

the std::mt19937 is not thread safe but you use the global mt from each thread without synchronization. Instead declare an initialize it in threadOptimize to give each thread its own RNG.

Same with bestSequence surround each copy to and from it with a mutex guard like you do when printing out.

The following for loop doesn't quite do what you think it does.

for(int i = 0; i  < distDim(mt);++i)
    seq[distDim(mt)] ^= true;

the distDim in the condition is reevaluated each time so it will likely never reach the max amount of iterations. You can solve that by going in reverse order:

for(int i = distDim(mt); i  > 0;--i)
    seq[distDim(mt)] ^= true;

But even then it won't be a perfect manipulation.

\$\endgroup\$
  • \$\begingroup\$ I edited my program with your suggestions. For that loop part I don't really care if it doesn't do what i think it does. .) i just wanted random behaviour. \$\endgroup\$ – TomazStoiljkovic Dec 3 '15 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.