for (auto j = i + 1; j < ds.size(); ++j) {

            vals[i][j] = 1.0 - std::exp(-(ds[i] * ds[j] * factor));

        }

    }

    double result = 0.0;

    for (auto i = 0u; i < ds.size(); ++i) {

        for (auto j = i + 1; j < ds.size(); ++j) {

            for (auto k = j + 1; k < ds.size(); ++k) {

                result += vals[i][j] * vals[j][k] * vals[i][k];

            }

        }

    }

    return result;

}

int main(int argc, char* argv[]) {

    // Simulation parameters

    const int numVerticesMin = 100;

    const int numVerticesMax = 1000;

    const int numVerticesStep = 100;

    float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f};

    auto getMixingTime = [](int n, float tau) {

    };

    constexpr int measurements = 10;

    constexpr int measureSkip = 1000; // Take a sample every ... steps

    // Output file

    std::ofstream outfile;

    if (argc >= 2)

        outfile.open(argv[1]);

    else

        outfile.open("graphdata_properties.m");

    if (!outfile.is_open()) {

        std::cout << "ERROR: Could not open output file.\n";

        return 1;

    }

    // Output Mathematica-style comment to indicate file contents

    outfile << "(*\n";

    outfile << "n from " << numVerticesMin << " to " << numVerticesMax

            << " step " << numVerticesStep << std::endl;

    outfile << "tauValues: " << tauValues << std::endl;

    outfile << "degreeSamples: " << totalDegreeSamples << std::endl;

    outfile << "data:\n";

    outfile << "1: {n,tau}\n";

    outfile << "2: avgTriangles\n";

    outfile << "3: edges\n";

    outfile << "4: dstn\n";

    outfile << "5: { HH A, HH L, average A, average L }  where for each there is (average of) {lambda1 , lambda1 - lambda2, lambda1/lambda2}\n";

    outfile << "6: switching successrate after mixing\n";

    outfile << "7: initial HH triangles\n";

    outfile << "*)" << std::endl;

    // Mathematica does not accept normal scientific notation

    outfile << std::fixed;

    outfile << '{';

    bool outputComma = false;

    std::mt19937 rng(std::random_device{}());

    Graph g;

    for (int numVertices = numVerticesMin; numVertices <= numVerticesMax;

         numVertices += numVerticesStep) {

        for (float tau : tauValues) {

            // For a single n,tau take samples over several instances of

            // the degree distribution.

            for (int degreeSample = 0; degreeSample < totalDegreeSamples;

                 ++degreeSample) {