Initial commit

Compile the code with an intel compiler (with MKL available):

mpiicpc -O3 -std=c++11 -mkl -I./ main.cpp -o rfm_code

The code itself can be run by

./rfm_code steps=4 type=h weik=2 niter=10000

which starts a network of size 2**4, network of type hierarchical

(nonhierarchical would be s) up to 10000 iterations (or failure). Niter

should be set large enough to be longer than the simulation duration. 

Running the simulations is in run.py and the statistical information needed for the creep simulations are calculated by collect_data.py. Figure S1 can be created by fig-s1.py, but it the statistical information needs to be calculated first for that.

If you want to make sure that it works (in the sense that the program does not fail/break down) on your system, execute test.py. If you want to retrace the paper, execute redo, but be aware that it might take weeks/months/year on a simple work station, as we used a high perofmrance cluster for that.  

#ifndef __RANDOM_H__

#define __RANDOM_H__

typedef unsigned long int U_Int_32_t;

typedef unsigned long long int  U_Int_64_t;

class Random {

public:

    U_Int_64_t u,v,w;

    Random(U_Int_64_t j) : v(4101842887655102017LL), w(1) {

        u = j ^ v; int64();

        v = u; int64();

        w = v; int64();

    }

    inline U_Int_64_t int64() {

        u = u * 2862933555777941757LL + 7046029254386353087LL;

        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;

        w = 4294957665U*(w & 0xffffffff) + (w >> 32);

        U_Int_64_t x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;

        return (x + v) ^ w;

    }

    inline double doub() { return 5.42101086242752217E-20 * int64(); }

    inline unsigned int int32() { return (unsigned int)int64(); }

};

#endif

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Convert network of diamond shape to image in order to be able

to use image convolutions etc. on it.

@author: Stefan Hiemer

"""

import os

import glob

from shutil import copy2

import numpy as np

import subprocess

from multiprocessing import cpu_count

from functools import partial

import json

from itertools import product

def check_dir(dir_name):

    """

    Checks wether directory exists. If true, moves in, else creates it and

    moves in.

    """

    if os.path.exists(dir_name):

        os.chdir(dir_name)

    else:

        os.mkdir(dir_name)

        os.chdir(dir_name)

    return

def get_data(sizes, ks, network_types):

    """

    Collect peak current and voltage at peak current from quasistatic random

    fuse simulations.

    sizes: list of ints

    ks: list of floats

    network_types: str, either 's' or 'h' which stands for nonhierarchical and

                   hierarchical.

    """

    home = os.getcwd()

    nproc = cpu_count() 

    # set up dictionary to collect the data

    data = dict()

    for network_type in network_types:

        data[network_type] = dict()

        for size in sizes:

            data[network_type]["step-"+str(size)] = dict()

            for k in ks:

                data[network_type]["step-"+str(size)]["k-"+str(k)] = dict()

                data[network_type]["step-"+str(size)]["k-"+str(k)]["seeds"] = []

                data[network_type]["step-"+str(size)]["k-"+str(k)]["imax"] = []

                data[network_type]["step-"+str(size)]["k-"+str(k)]["vatmax"] = []

    for network_type,size,k in product(network_types,sizes,ks):

        if network_type == 'h':

            os.chdir("hierarchical")

        elif network_type == 's':

            os.chdir("nonhierarchical")

        else:

            raise ValueError("Unknown network type")

        os.chdir("step-"+str(size))

        os.chdir("k-"+str(k))

        # get seeds

        seeds = glob.glob("seed-*")

        i = 0

        for seed in seeds:

            os.chdir(seed)

            # get number of iv_seeds

            files = glob.glob("iv_seed*")

            for file in files:

                print(i+1,network_type, size, k, seed, file)

                i+=1

                # load file

                vi = np.loadtxt(file,dtype=np.float,delimiter=' ')

                #save seed

                data[network_type]["step-"+str(size)]["k-"+str(k)]["seeds"].append(int(''.join(file.rsplit("_",1))[7:]))

                # get maximum threshold

                ind = np.argmax(vi[:,1])

                data[network_type]["step-"+str(size)]["k-"+str(k)]["vatmax"].append(vi[ind,0])

                data[network_type]["step-"+str(size)]["k-"+str(k)]["imax"].append(vi[ind,1])

            os.chdir("../")

        # move out of folders after job is finished

        os.chdir("../../../")

    os.chdir(home)

    with open("data.json", "w") as file:

        json.dump(data,file)

    return

def stat_moments(sizes, ks, network_types):

    """

    Calculate mean, median and standard deviation for maximum/peak current and

    voltage at peak current.

    sizes: list of ints

    ks: list of floats

    network_types: str, either 's' or 'h' which stands for nonhierarchical and

                   hierarchical.

    """

    # set up dictionary to collect the data

    scale = dict()

    for network_type in network_types:

        scale[network_type] = dict()

        for size in sizes:

            scale[network_type]["step-"+str(size)] = dict()

            for k in ks:

                scale[network_type]["step-"+str(size)]["k-"+str(k)] = dict()

    # load data from json

    with open("data.json", "r") as file:

        data = json.load(file)

    for network_type,size,k in product(network_types,sizes,ks):

        imax = data[network_type]["step-"+str(size)]["k-"+str(k)]["imax"]

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["av_imax"] = np.mean(imax)

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["med_imax"] = np.median(imax)

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["std_imax"] = np.std(imax)

        vatmax = data[network_type]["step-"+str(size)]["k-"+str(k)]["vatmax"]

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["av_vatmax"] = np.mean(vatmax)

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["med_vatmax"] = np.median(vatmax)

        scale[network_type]["step-"+str(size)]["k-"+str(k)]["std_vatmax"] = np.std(vatmax)

    with open("stat-moments.json", "w") as file:

        json.dump(scale,file)

    return

def convergence(splits,sizes, ks, network_types):

    # set up dictionary to collect the data

    deviation = dict()

    for network_type in network_types:

        deviation[network_type] = dict()

        for size in sizes:

            deviation[network_type]["step-"+str(size)] = dict()

            for k in ks:

                deviation[network_type]["step-"+str(size)]["k-"+str(k)] = dict()

                deviation[network_type]["step-"+str(size)]["k-"+str(k)]["samples"] = []

                deviation[network_type]["step-"+str(size)]["k-"+str(k)]["std"] = []

                deviation[network_type]["step-"+str(size)]["k-"+str(k)]["var_coeff"] = []

    # load data from json

    with open("data.json", "r") as file:

        data = json.load(file)

    for network_type,size,k in product(network_types,sizes,ks):

        mean = np.mean(data[network_type]["step-"+str(size)]["k-"+str(k)]["imax"])

        for split in splits:

            reshaped = np.array(data[network_type]["step-"+str(size)]["k-"+str(k)]["imax"])

            reshaped = np.reshape(reshaped,

                                  (split,int(reshaped.shape[0]/split)))

            n = np.shape(reshaped)[1]

            std = np.std(np.mean(reshaped,axis=1))

            deviation[network_type]["step-"+str(size)]["k-"+str(k)]["samples"].append(n)

            deviation[network_type]["step-"+str(size)]["k-"+str(k)]["std"].append(std)

            deviation[network_type]["step-"+str(size)]["k-"+str(k)]["var_coeff"].append(std/mean)

    with open("deviation.json", "w") as file:

        json.dump(deviation,file)

    return

if __name__ == "__main__":

    get_data(sizes=[9,8,7,6,5,4,3,2],

             ks=[1,2,4,8,16],

             network_types=['s','h'])

    stat_moments(sizes=[9,8,7,6,5,4,3,2],

                 ks=[1,2,4,8,16],

                 network_types=['s','h'])

#ifndef EXECUTIONBLOCK_HPP

#define EXECUTIONBLOCK_HPP

#include <sstream>

#include <iostream>

#include <outiv.hpp>

#include <Random.h>

//#include <lattice.hpp>

#include <lattice_shufflings.hpp>

#include <fuser.hpp>

class ExecutionBlock {

    uint64_t seed;

    int steps;

    char type;

    double weik;

    int nremove;

    int niter;

    double Ifixed;

    double temp;

    int networks;

    bool optional;

    bool storeiv;

    public:

    ExecutionBlock(

            uint64_t seed,

            int steps,

            char type,

            double weik,

            int nremove,

            int niter,

            //double Ifixed,

            //double temp,

            int networks,

            bool optional,

            bool storeiv) :

        seed(seed),

        steps(steps),

        type(type),

        weik(weik),

        nremove(nremove),

        niter(niter),

        //Ifixed(Ifixed),

        //temp(temp),

        networks(networks),

        optional(optional),

        storeiv(storeiv) {}

    OutIV run() {

        Random ra(seed);

        int realizations = networks;

        OutIV outIV(realizations);

        //if(realizations>1) optional = false;

        for(int r=0; r<realizations; ++r) {

            Lattice l(ra,steps,type,weik,nremove);

            if(optional) {

                stringstream e_stream;

                e_stream << "edgelist_0_seed" << seed << "_" << r;

                l.export_edgelist_nobus(e_stream.str().c_str());

                //l.export_nodes_coo("nodes");

                //l.export_links_coo("links");

                //l.export_links_extremes("segments");

            }

            Fuser f(l/*,Ifixed,temp*/);

            //f.execute(niter);

            //outIV.accumulate_envelope(f.execute(niter), optional);

            //outIV.accumulate_avalanches(f.execute(niter));

            f.execute(niter); //this choice circumvents outIV completely

            if(storeiv) {

                stringstream namestream;

                namestream << "iv_seed" << seed << "_" << r;

                //f.export_qstatic_iv_only(namestream.str().c_str());

                f.export_qstatic_iv_only_doubleprecision(namestream.str().c_str());

                //namestream << "tvp_seed" << seed << "_" << r;

                //f.export_qstatic_tVp_singleprecision(namestream.str().c_str());

                //f.export_qstatic_tVp_doubleprecision(namestream.str().c_str());

                //stringstream avastream;

                //avastream << "avasizes_seed" << seed << "_" << r;

                //f.export_qstatic_avasizes_singleprecision(avastream.str().c_str());

            }

            stringstream p_stream;

            p_stream << "potentials_seed" << seed << "_" << r;

            f.export_potentials_nodewise(p_stream.str().c_str());

            //stringstream b_stream;

            //b_stream << "brokenlist_seed" << seed << "_" << r;

            //l.export_brokenlist_nobus(b_stream.str().c_str());

            stringstream bs_stream;

            bs_stream << "brokenlist_sorted_seed" << seed << "_" << r;

            l.export_brokenlist_sorted_nobus(bs_stream.str().c_str());

            if(optional) {

                //l.export_edgelist_nobus("edgelist_final");

                //l.export_broken_coo("broken");

                //l.export_links_coo("surviving");

                //f.export_potentials("potentials");

                //f.export_vcurrents("vcurrents");

                //stringstream namestream;

                //namestream << "crack_seed" << seed << "_" << r;

                //f.export_crack(namestream.str().c_str());

                //l.export_links_extremes("surviving");

                //f.export_qstatic_iv("iv");

                //cout << "pbroken " << l.nbroken/(double)l.N << endl;

            }

        }

        //outIV.export_average_envelope();

        return outIV;

    }

};

#endif