TARGETS += switchchain_canonical_properties

TARGETS += switchchain_exponent_new

void generateCanonicalPowerlawGraph(int n, float tau, Graph &g,

                                    DegreeSequence &ds) {

    ds.resize(n);

    powerlaw_distribution degDist(tau, 1, n);

    degDist.getFixedDistribution(n, ds.begin());

    unsigned int sum = std::accumulate(ds.begin(), ds.end(), 0);

    if (sum % 2) {

        //TODO: Can we do this: ??

        ds.back()++;

    }

    if (g.createFromDegreeSequence(ds))

        return;

    g.reset(n);

    std::cerr << "ERROR: Canonical ds is not graphical" << std::endl;

}

#include <iostream>

// Let m = minimum and M = maximum and a = alpha

// Then for continuous version we have

// PDF: f(x) = C x^{-a}

// with C = (1-a)/(M^{1-a} - m^{1-a})

// This gives CDF: F(x) = (x^{1-a} - m^{1-a}) / (M^{1-a} - m^{1-a})

// with inversion F^{-1}(y) = ( y M^{1-a} + (1-y) m^{1-a} )^{1/(1-a)}

// i.e. linear interpolate between M^{1-a} < m^{1-a}

// For M = infty this becomes

//  F'^{-1}(y) = m (1-y)^{1/(1-a)}

//             = m (y')^{-1/(a-1)}

// where higher y' means lower x

    // Generate non-random ``canonical'' powerlaw distribution

    // Same as above operator but where the random number between

    // 0 and 1 is replaced by stepping from 0 to 1 in fixed steps

    template <class FwdIterator>

    void getFixedDistribution(int n, FwdIterator output) const {

        // go in linear steps over interval

        // [ M^{1-a} , m^{1-a} ]

        double minVal = std::pow(double(xmax), 1.0f - alpha);

        double maxVal = std::pow(double(xmin), 1.0f - alpha);

        unsigned int x;

        // Now it outputs in increasing order

        for (int i = 0; i < n; ++i) {

            double r1 = double(i) / double(n);

            double r2 = r1 * minVal + (1.0 - r1) * maxVal;

            x = std::round(std::pow(r2, _exponent));

            if (x > xmax || x < xmin) {

                std::cerr << "ERROR: x not in [xmin,xmax] : " << x

                          << " not in [" << xmin << ',' << xmax

                          << "]; i = " << i << std::endl;

            }

            *output++ = x;

        }

    }

#include "exports.hpp"

#include "graph.hpp"

#include "graph_powerlaw.hpp"

#include "graph_spectrum.hpp"

#include "switchchain.hpp"

#include <algorithm>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

double getDSTN(const DegreeSequence& ds) {

    std::vector<std::vector<double>> vals(ds.size());

    for (auto& v : vals) {

        v.resize(ds.size(), 0);

    }

    auto D = 0u;

    for (auto d : ds)

        D += d;

    double factor = 1.0 / double(D);

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

        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 = 1000;

    const int numVerticesMax = 10000;

    const int numVerticesStep = 1000;

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

    //const int totalDegreeSamples = 5000;

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

        return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n);

    };

    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_canonical_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 << "Canonical degree sequence.\n";

    outfile << "mixingTime: 50 * (50 - 30 (tau - 2)) n\n";

    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;

    Graph g;

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

         numVertices += numVerticesStep) {

        for (float tau : tauValues) {

            DegreeSequence ds;

            generateCanonicalPowerlawGraph(numVertices, tau, g, ds);

            SwitchChain chain;

            if (!chain.initialize(g)) {

                std::cerr << "Could not initialize Markov chain.\n";

                return 1;

            }

            std::cout << "Running (n,tau) = (" << numVertices << ',' << tau

                      << "). " << std::flush;

            // Mix

            int mixingTime = getMixingTime(numVertices, tau);

            for (int i = 0; i < mixingTime; ++i) {

                chain.doMove();

            }

            std::cout << "Mixing done. " << std::flush;

            std::array<double, 3> HHAspectrum;

            std::array<double, 3> HHLspectrum;

            std::array<double, 3> avgAspectrum;

            std::array<double, 3> avgLspectrum;

            auto getSpectralValues =

                [](const std::vector<float> &s) -> std::array<double, 3> {

                auto l1 = s[s.size() - 1];

                auto l2 = s[s.size() - 2];

                return {l1, l1 - l2, l1 / l2};

            };

            GraphSpectrum gs_start(g);

            GraphSpectrum gs(chain.g);

            HHAspectrum =

                getSpectralValues(gs_start.computeAdjacencySpectrum());

            HHLspectrum =

                getSpectralValues(gs_start.computeLaplacianSpectrum());

            long long trianglesTotal = 0;

            int movesDone = 0;

            avgAspectrum.fill(0);

            avgLspectrum.fill(0);

            for (int i = 0; i < measurements; ++i) {

                for (int j = 0; j < measureSkip; ++j)

                    if (chain.doMove())

                        ++movesDone;

                trianglesTotal += chain.g.countTriangles();

                auto sA = getSpectralValues(gs.computeAdjacencySpectrum());

                auto sL = getSpectralValues(gs.computeLaplacianSpectrum());

                for (auto i = 0u; i < 3; ++i) {

                    avgAspectrum[i] += sA[i];

                    avgLspectrum[i] += sL[i];

                }

            }

            float avgTriangles = float(trianglesTotal) / float(measurements);

            float successrate =

                float(movesDone) / float(measurements * measureSkip);

            for (auto &f : avgAspectrum)

                f /= float(measurements);

            for (auto &f : avgLspectrum)

                f /= float(measurements);

            std::cout << "Measuring done." << std::flush;

            if (outputComma)

                outfile << ',' << '\n';

            outputComma = true;

            outfile << '{' << '{' << numVertices << ',' << tau << '}';

            outfile << ',' << avgTriangles;

            outfile << ',' << g.edgeCount();

            outfile << ',' << getDSTN(ds);

            outfile << ',' << '{' << HHAspectrum;

            outfile << ',' << HHLspectrum;

            outfile << ',' << avgAspectrum;

            outfile << ',' << avgLspectrum;

            outfile << '}';

            outfile << ',' << successrate;

            outfile << ',' << g.countTriangles();

            outfile << '}' << std::flush;

            std::cout << "Output done." << std::endl;

        }

    }

    outfile << '}';

    return 0;

}

#include "exports.hpp"

#include "graph.hpp"

#include "graph_powerlaw.hpp"

#include "graph_spectrum.hpp"

#include "switchchain.hpp"

#include <algorithm>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

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

    // Simulation parameters

    const int numVerticesMin = 20000;

    const int numVerticesMax = 50000;

    const int numVerticesStep = 10000;

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

    float tauValues[] = {2.1f, 2.3f, 2.5f, 2.7f, 2.9f};

    const int totalDegreeSamples = 1000;

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

        return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n);

    };

    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_exponent_new.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 << "mixingTime: 50 * (50 - 30 (tau - 2)) n\n";

    outfile << "data:\n";

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

    outfile << "2: avgTriangles\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) {

                DegreeSequence ds;

                generatePowerlawGraph(numVertices, tau, g, ds, rng);

                SwitchChain chain;

                if (!chain.initialize(g)) {

                    std::cerr << "Could not initialize Markov chain.\n";

                    return 1;

                }

                std::cout << "Running (n,tau) = (" << numVertices << ','

                          << tau << "). " << std::flush;

                // Mix

                int mixingTime = getMixingTime(numVertices, tau);

                for (int i = 0; i < mixingTime; ++i) {

                    chain.doMove();

                }

                std::cout << "Mixing done. " << std::flush;

                long long trianglesTotal = 0;

                for (int i = 0; i < measurements; ++i) {

                    for (int j = 0; j < measureSkip; ++j)

                        chain.doMove();

                    trianglesTotal += chain.g.countTriangles();

                }

                float avgTriangles =

                    float(trianglesTotal) / float(measurements);

                std::cout << "Measuring done." << std::flush;

                if (outputComma)

                    outfile << ',' << '\n';

                outputComma = true;

                outfile << '{' << '{' << numVertices << ',' << tau << '}';

                outfile << ',' << avgTriangles;

                outfile << '}' << std::flush;

                std::cout << "Output done." << std::endl;

            }

        }

    }

    outfile << '}';

    return 0;

}

(* When importing from exponent-only-data file or property-data file *)

(* graphdata_exponent_mix32.m *)

(* graphdata_exponent_highN.m *)

(* graphdata_properties2.m *)

nRange=5;;-1;

nRange=All;

loglogdata=Log[averagesErrorBars[[All,nRange,1]]];

gsraw=Import[NotebookDirectory[]<>"data/graphdata_properties.m"];

2: avgTriangles

3: edges

4: dstn

5: { HH A, HH L, average A, average L }  where for each there is (average of) {lambda1 , lambda1 - lambda2, lambda1/lambda2}

6: switching successrate after mixing

7: initial HH triangles

Histogram[gdata[[1,1,All,5,3,1]]]

Histogram[gdata[[1,1,All,5,3,2]]]