push

README.md

 # ContinuousIM
 
+Code for experiments of Chapter 3 of the thesis 'Influence Maximisation and Adaptive Control of Opinion Dynamics on Complex Networks'.
 
-
-## Getting started
-
-To make it easy for you to get started with GitLab, here's a list of recommended next steps.
-
-Already a pro? Just edit this README.md and make it your own. Want to make it easy? [Use the template at the bottom](#editing-this-readme)!
-
-## Add your files
-
-- [ ] [Create](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#create-a-file) or [upload](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#upload-a-file) files
-- [ ] [Add files using the command line](https://docs.gitlab.com/ee/gitlab-basics/add-file.html#add-a-file-using-the-command-line) or push an existing Git repository with the following command:
-
-```
-cd existing_repo
-git remote add origin https://git.soton.ac.uk/grm1g17/continuousim.git
-git branch -M main
-git push -uf origin main
-```
-
-## Integrate with your tools
-
-- [ ] [Set up project integrations](https://git.soton.ac.uk/grm1g17/continuousim/-/settings/integrations)
-
-## Collaborate with your team
-
-- [ ] [Invite team members and collaborators](https://docs.gitlab.com/ee/user/project/members/)
-- [ ] [Create a new merge request](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html)
-- [ ] [Automatically close issues from merge requests](https://docs.gitlab.com/ee/user/project/issues/managing_issues.html#closing-issues-automatically)
-- [ ] [Enable merge request approvals](https://docs.gitlab.com/ee/user/project/merge_requests/approvals/)
-- [ ] [Set auto-merge](https://docs.gitlab.com/ee/user/project/merge_requests/merge_when_pipeline_succeeds.html)
-
-## Test and Deploy
-
-Use the built-in continuous integration in GitLab.
-
-- [ ] [Get started with GitLab CI/CD](https://docs.gitlab.com/ee/ci/quick_start/index.html)
-- [ ] [Analyze your code for known vulnerabilities with Static Application Security Testing(SAST)](https://docs.gitlab.com/ee/user/application_security/sast/)
-- [ ] [Deploy to Kubernetes, Amazon EC2, or Amazon ECS using Auto Deploy](https://docs.gitlab.com/ee/topics/autodevops/requirements.html)
-- [ ] [Use pull-based deployments for improved Kubernetes management](https://docs.gitlab.com/ee/user/clusters/agent/)
-- [ ] [Set up protected environments](https://docs.gitlab.com/ee/ci/environments/protected_environments.html)
-
-***
-
-# Editing this README
-
-When you're ready to make this README your own, just edit this file and use the handy template below (or feel free to structure it however you want - this is just a starting point!). Thank you to [makeareadme.com](https://www.makeareadme.com/) for this template.
-
-## Suggestions for a good README
-Every project is different, so consider which of these sections apply to yours. The sections used in the template are suggestions for most open source projects. Also keep in mind that while a README can be too long and detailed, too long is better than too short. If you think your README is too long, consider utilizing another form of documentation rather than cutting out information.
-
-## Name
-Choose a self-explaining name for your project.
-
-## Description
-Let people know what your project can do specifically. Provide context and add a link to any reference visitors might be unfamiliar with. A list of Features or a Background subsection can also be added here. If there are alternatives to your project, this is a good place to list differentiating factors.
-
-## Badges
-On some READMEs, you may see small images that convey metadata, such as whether or not all the tests are passing for the project. You can use Shields to add some to your README. Many services also have instructions for adding a badge.
-
-## Visuals
-Depending on what you are making, it can be a good idea to include screenshots or even a video (you'll frequently see GIFs rather than actual videos). Tools like ttygif can help, but check out Asciinema for a more sophisticated method.
-
-## Installation
-Within a particular ecosystem, there may be a common way of installing things, such as using Yarn, NuGet, or Homebrew. However, consider the possibility that whoever is reading your README is a novice and would like more guidance. Listing specific steps helps remove ambiguity and gets people to using your project as quickly as possible. If it only runs in a specific context like a particular programming language version or operating system or has dependencies that have to be installed manually, also add a Requirements subsection.
-
-## Usage
-Use examples liberally, and show the expected output if you can. It's helpful to have inline the smallest example of usage that you can demonstrate, while providing links to more sophisticated examples if they are too long to reasonably include in the README.
-
-## Support
-Tell people where they can go to for help. It can be any combination of an issue tracker, a chat room, an email address, etc.
-
-## Roadmap
-If you have ideas for releases in the future, it is a good idea to list them in the README.
-
-## Contributing
-State if you are open to contributions and what your requirements are for accepting them.
-
-For people who want to make changes to your project, it's helpful to have some documentation on how to get started. Perhaps there is a script that they should run or some environment variables that they need to set. Make these steps explicit. These instructions could also be useful to your future self.
-
-You can also document commands to lint the code or run tests. These steps help to ensure high code quality and reduce the likelihood that the changes inadvertently break something. Having instructions for running tests is especially helpful if it requires external setup, such as starting a Selenium server for testing in a browser.
-
-## Authors and acknowledgment
-Show your appreciation to those who have contributed to the project.
-
-## License
-For open source projects, say how it is licensed.
-
-## Project status
-If you have run out of energy or time for your project, put a note at the top of the README saying that development has slowed down or stopped completely. Someone may choose to fork your project or volunteer to step in as a maintainer or owner, allowing your project to keep going. You can also make an explicit request for maintainers.

datasets/classroom/Friendship-network_data_2013.csv

+1 55
+1 205
+1 272
+1 494
+1 779
+1 894
+3 1
+3 28
+3 147
+3 272
+3 407
+3 674
+3 884
+27 63
+27 173
+28 202
+28 327
+28 353
+28 407
+28 429
+28 441
+28 492
+28 545
+32 440
+32 624
+32 797
+32 920
+34 151
+34 277
+34 502
+34 866
+45 48
+45 79
+45 335
+45 496
+45 601
+45 674
+45 765
+46 117
+46 196
+46 257
+46 268
+48 45
+48 79
+48 496
+55 1
+55 170
+55 205
+55 252
+55 272
+55 779
+55 883
+55 894
+61 797
+63 27
+63 125
+63 173
+70 101
+70 132
+70 240
+70 425
+70 447
+72 407
+72 674
+72 857
+79 45
+79 48
+79 335
+79 496
+79 601
+79 674
+79 765
+80 120
+80 285
+80 468
+80 601
+85 190
+85 213
+85 214
+85 335
+85 603
+85 605
+85 765
+92 468
+92 845
+101 70
+101 119
+101 122
+101 132
+101 240
+101 343
+101 364
+101 425
+101 447
+117 1
+117 46
+117 196
+117 205
+117 252
+117 257
+117 265
+117 268
+117 272
+117 364
+117 407
+117 465
+117 494
+117 587
+117 883
+117 894
+119 101
+119 122
+119 132
+119 240
+119 425
+119 447
+120 80
+120 285
+120 488
+122 101
+122 119
+122 255
+122 425
+122 447
+124 471
+124 970
+125 92
+125 248
+125 325
+125 468
+125 491
+125 622
+125 624
+125 797
+125 960
+132 70
+132 101
+132 119
+132 122
+132 240
+132 425
+132 447
+134 388
+134 492
+134 496
+147 1
+147 3
+147 28
+147 72
+147 184
+147 272
+147 353
+147 407
+147 674
+147 884
+151 34
+151 38
+151 201
+151 277
+151 452
+151 502
+151 634
+151 642
+151 691
+151 694
+151 753
+151 866
+151 869
+156 694
+159 38
+159 642
+165 498
+170 1
+170 55
+170 205
+170 779
+170 883
+170 894
+173 27
+173 63
+173 1332
+184 327
+184 429
+184 441
+190 85
+190 213
+190 214
+190 272
+196 46
+196 117
+196 252
+196 268
+196 364
+200 480
+200 845
+201 34
+201 38
+201 245
+201 502
+201 642
+201 691
+201 753
+201 869
+202 28
+202 170
+202 407
+202 545
+202 883
+205 1
+205 55
+205 117
+205 170
+205 265
+205 272
+205 494
+205 587
+205 691
+205 779
+205 883
+205 894
+211 242
+211 468
+211 845
+213 190
+213 214
+214 190
+214 213
+214 603
+219 605
+222 248
+222 343
+222 867
+232 492
+232 798
+240 70
+240 101
+240 119
+240 132
+240 327
+240 425
+240 447
+242 211
+242 468
+242 845
+245 325
+245 440
+245 634
+245 691
+245 753
+245 869
+245 959
+248 222
+248 564
+248 694
+252 55
+252 117
+252 196
+252 265
+252 272
+252 364
+252 779
+255 275
+257 46
+257 268
+257 364
+265 117
+265 170
+265 196
+265 205
+265 252
+265 494
+265 587
+265 883
+265 894
+268 46
+268 117
+268 196
+268 257
+268 364
+272 1
+272 55
+272 170
+272 190
+272 205
+272 214
+272 441
+272 494
+272 587
+272 779
+272 883
+275 255
+275 312
+275 612
+277 34
+277 151
+277 502
+277 634
+277 691
+277 866
+285 80
+285 120
+285 232
+285 488
+285 492
+312 275
+312 612
+325 125
+325 245
+325 622
+325 624
+325 769
+325 797
+325 959
+327 27
+327 28
+327 119
+327 184
+327 353
+327 407
+327 429
+327 441
+335 45
+335 79
+335 765
+343 101
+343 222
+343 867
+353 28
+353 46
+353 122
+353 407
+353 425
+364 117
+364 196
+364 252
+364 257
+364 268
+366 974
+388 45
+388 79
+388 134
+388 335
+388 492
+388 496
+388 603
+388 765
+407 3
+407 28
+407 72
+407 147
+407 184
+407 202
+407 327
+407 353
+407 429
+407 441
+407 545
+407 674
+407 884
+425 70
+425 101
+425 119
+425 122
+425 132
+425 240
+425 343
+425 353
+425 441
+425 447
+429 28
+429 119
+429 184
+429 327
+429 353
+429 407
+429 441
+440 32
+440 245
+440 605
+440 797
+440 920
+441 28
+441 184
+441 272
+441 327
+441 407
+441 429
+441 447
+447 70
+447 101
+447 119
+447 122
+447 132
+447 240
+447 425
+447 441
+452 85
+452 634
+452 691
+452 869
+452 1332
+465 486
+465 531
+465 857
+468 80
+468 92
+468 125
+468 211
+468 242
+468 601
+468 845
+471 124
+471 970
+480 200
+480 771
+486 465
+486 531
+488 48
+488 120
+488 285
+491 219
+491 520
+491 576
+491 605
+492 28
+492 120
+492 134
+492 232
+492 285
+492 388
+492 447
+492 488
+492 496
+494 1
+494 117
+494 205
+494 265
+494 272
+494 587
+494 883
+494 894
+496 45
+496 48
+496 79
+496 85
+496 134
+496 388
+496 492
+496 603
+498 165
+498 857
+502 151
+502 277
+502 691
+502 866
+502 869
+520 219
+520 491
+520 576
+520 605
+531 465
+531 486
+531 691
+545 28
+545 202
+545 407
+564 577
+564 694
+576 219
+576 491
+576 520
+576 605
+577 564
+577 694
+587 117
+587 265
+587 272
+587 494
+587 883
+601 45
+601 79
+601 80
+601 134
+601 200
+601 388
+601 603
+603 85
+603 214
+603 335
+603 388
+603 496
+603 765
+605 85
+605 201
+605 219
+605 520
+605 576
+612 275
+612 312
+622 125
+622 452
+622 624
+622 769
+622 797
+622 798
+622 959
+622 960
+624 125
+624 245
+624 325
+624 452
+624 491
+624 622
+624 769
+624 797
+624 798
+624 960
+634 691
+634 869
+634 1332
+642 151
+642 201
+642 452
+642 634
+642 691
+642 866
+642 869
+674 3
+674 45
+674 72
+674 79
+674 147
+674 407
+691 205
+691 245
+691 452
+691 502
+691 531
+691 634
+691 642
+691 869
+691 883
+691 1332
+694 564
+694 577
+753 201
+753 245
+765 45
+765 79
+765 335
+769 125
+769 245
+769 325
+769 452
+769 622
+769 624
+769 797
+769 959
+769 960
+771 200
+771 480
+779 1
+779 48
+779 55
+779 170
+779 205
+779 252
+779 272
+779 883
+779 894
+797 125
+797 325
+797 622
+797 624
+797 959
+798 605
+798 960
+845 80
+845 92
+845 200
+845 211
+845 242
+845 468
+857 55
+857 72
+857 465
+857 498
+857 779
+866 34
+866 151
+866 452
+866 502
+866 634
+866 642
+866 691
+866 753
+866 869
+867 222
+867 343
+869 245
+869 634
+869 691
+869 1332
+883 1
+883 55
+883 117
+883 170
+883 205
+883 252
+883 265
+883 272
+883 494
+883 587
+883 779
+883 894
+883 1401
+884 3
+884 272
+894 1
+894 55
+894 117
+894 170
+894 205
+894 265
+894 272
+894 494
+894 587
+894 779
+894 883
+920 32
+920 440
+920 797
+959 245
+959 325
+959 622
+959 624
+959 769
+959 797
+960 622
+960 798
+970 124
+970 471
+974 366
+974 1485
+1228 642
+1228 1401
+1228 1519
+1332 452
+1332 634
+1332 691
+1332 869
+1332 1401
+1332 1519
+1401 642
+1401 1228
+1401 1332
+1401 1519
+1485 974
+1519 642
+1519 1228
+1519 1332
+1519 1401
+1519 1594
+1519 1828
+1594 1519
+1594 1828
+1828 1519
+1828 1594

datasets/classroom/output.out

+\BOOKMARK [1][-]{section.1}{Details of the network employed}{}% 1
+\BOOKMARK [1][-]{section.2}{First-order Taylor expansion of the vote share for small allocations}{}% 2
+\BOOKMARK [1][-]{section.3}{Zero-order Taylor expansion of the vote share for large allocations}{}% 3
+\BOOKMARK [1][-]{section.4}{Exploitability of a passive controller when deviating from uniform allocation}{}% 4
+\BOOKMARK [1][-]{section.5}{Comparison of numerical and analytical results for optimal allocations in the continuous regime and budget disadvantage}{}% 5
+\BOOKMARK [1][-]{section.6}{Shadowing, shielding and hub preferences in the discrete regime}{}% 6
+\BOOKMARK [1][-]{section.7}{Proof of uniqueness of Nash Equilibrium}{}% 7
+\BOOKMARK [1][-]{section.8}{Numerical results of an iterative approach to finding the Nash Equilibrium}{}% 8
+\BOOKMARK [1][-]{section.9}{Extension to other network topologies}{}% 9
+\BOOKMARK [1][-]{section.10}{Proof of concavity}{}% 10
+\BOOKMARK [1][-]{section.11}{Iterations to optimal IM of the proposed heuristics}{}% 11

datasets/ia-email-univ/readme.html

+<html>
+<link rel="stylesheet" media="screen" href="http://networkrepository.com/assets/css/bootstrap.css">
+<link rel="stylesheet" href="http://networkrepository.com/assets/css/responsive.css">   
+<link rel="stylesheet" href="http://networkrepository.com/assets/css/style.css">
+
+<body>
+<br><br>
+<div class="container">
+<div class="row">
+<div class="col-md-12">
+<div class="panel panel-success">
+    <div class="panel-heading">
+        <h3 class="panel-title"><i class="icon-group icon-large"></i> Acknowledgement policy</h3>
+    </div>
+    <div class="panel-body">
+        <h4>Please acknowledge the <a href="http://networkrepository.com" target="_blank">repository</a> in published materials</h4> 
+        <p class="lead12">                    
+            If you publish material based on data obtained from this <a href="http://networkrepository.com" target="_blank">repository</a>, then, in your acknowledgements, please note the assistance you received by using <a href="http://networkrepository.com" target="_blank">network repository</a>. 
+            <BR><BR>
+            Please use the following <code>BiBTeX</code> reference:<BR>
+
+            <blockquote class="hero">
+            <p style="font-family:'Lato';">
+            <code style="font-size:1.1em;">
+                        &nbsp; @inproceedings{nr-aaai15,<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    title = {The Network Data Repository with Interactive Graph Analytics and Visualization},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    author={Ryan A. Rossi and Nesreen K. Ahmed},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    booktitle = {Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    url={<i><a href="http://networkrepository.com" target="_blank"><b style="font-weight:bold !important;">http://networkrepository.com</b></a></i>},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    year={2015}<BR>
+                        &nbsp; }<BR>
+            </code>
+            </p>
+            </blockquote>
+            Many of the datasets have additional citation requests; these can be found on each dataset page.
+            See the <a href="http://networkrepository.com/policy.php" target="_blank">data license and policy</a> for more information.
+        </p>
+    </div>
+</div>
+</div>
+</div>
+</div>
+</body>
+</html>
\ No newline at end of file

datasets/soc-advogato/readme.html

+<html>
+<link rel="stylesheet" media="screen" href="http://networkrepository.com/assets/css/bootstrap.css">
+<link rel="stylesheet" href="http://networkrepository.com/assets/css/responsive.css">   
+<link rel="stylesheet" href="http://networkrepository.com/assets/css/style.css">
+
+<body>
+<br><br>
+<div class="container">
+<div class="row">
+<div class="col-md-12">
+<div class="panel panel-success">
+    <div class="panel-heading">
+        <h3 class="panel-title"><i class="icon-group icon-large"></i> Acknowledgement policy</h3>
+    </div>
+    <div class="panel-body">
+        <h4>Please cite the following if you use the data:</h4>
+        <p class="lead12">           
+
+            <blockquote class="hero">
+            <p style="font-family:'Lato';">
+            <code style="font-size:1.1em;">
+                        &nbsp; @inproceedings{nr-aaai15,<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    title = {The Network Data Repository with Interactive Graph Analytics and Visualization},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    author={Ryan A. Rossi and Nesreen K. Ahmed},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    booktitle = {AAAI},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    url={<i><a href="http://networkrepository.com" target="_blank"><b style="font-weight:bold !important;">http://networkrepository.com</b></a></i>},<BR>
+                        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;    year={2015}<BR>
+                        &nbsp; }<BR>
+            </code>
+            </p>
+            </blockquote>
+        </p>
+    </div>
+</div>
+</div>
+</div>
+</div>
+</body>
+</html>
\ No newline at end of file

hillclimb.py

+import time
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+PLOTS_PATH = 'plots/'
+
+
+class HillClimber:
+    def __init__(self, W, b_vec, n_nodes, gain, initial_a_disc="random", **pars):
+        self.W = W
+        self.N = W.shape[0]
+        self.L = np.diagflat(np.sum(self.W, axis=0)) - self.W
+        self.seq = np.squeeze(np.array(np.sum(self.W, axis=1)))
+
+        self.gain = gain
+        self.n_nodes = n_nodes
+        if self.n_nodes > self.N:
+            raise Exception(
+                "There are more nodes for discrete allocation ({:d}) than nodes in the network ({:d}).".format(
+                    self.n_nodes, self.N))
+        self.best_degrees = []  # For the greedy initialization
+
+        # Initialize allocations
+        self.b = b_vec
+        if type(initial_a_disc) == np.ndarray:
+            if len(np.unique(initial_a_disc)) == 2:
+                assert sum(initial_a_disc > 0) == n_nodes
+                self.a_ind = set(np.arange(self.N)[initial_a_disc > 0])
+                print("Initial 'a' from previous 'a'")
+            else:
+                assert np.isclose(np.sum(initial_a_disc), n_nodes * gain)
+                self.a_ind = set(np.argsort(initial_a_disc)[:n_nodes])
+                print("Initial 'a' from continuous 'a'")
+        elif initial_a_disc == "random":
+            self.a_ind = set(np.random.choice(self.N, self.n_nodes, replace=False))
+        elif initial_a_disc == "zeros":
+            X = self.greedy()
+            print("Hill climbing initialized with greedy algorithm (X={:f})".format(X))
+        else:
+            raise Exception("Initial a-allocation '" + str(initial_a_disc) + "' not understood.")
+
+    def produce_a_vec(self):
+        # self.a = np.zeros((self.N, 1))
+        self.a = np.zeros(self.N)
+        self.a[list(self.a_ind)] = self.gain
+        self.a = self.a[:, None]
+
+    def greedy(self):
+        self.a_ind = set()
+        self.produce_a_vec()
+
+        degree_dic = {}
+        for j, degree in enumerate(self.seq):
+            if degree not in degree_dic:
+                degree_dic[degree] = []
+            degree_dic[degree].append(j)
+
+        for k in range(self.n_nodes):
+            best_i, best_X, best_degree = None, -1, None
+            # for i, a_i in enumerate(self.a):
+            #     if np.isclose(a_i, 0):
+            #         self.a[i] = self.gain
+            #         X, x = self.compute_x()
+            #         if X > best_X:
+            #             best_X = X
+            #             best_i = i
+            #         self.a[i] = 0
+            for degree, indices in degree_dic.items():
+                i = np.random.choice(indices)
+                self.a[i] = self.gain
+                X, x = self.compute_x(produce_a_vec=False)
+                if X > best_X:
+                    best_X = X
+                    best_i = i
+                    best_degree = degree
+                self.a[i] = 0
+
+            self.a[best_i] = self.gain
+            self.a_ind.add(best_i)
+            degree_dic[best_degree].remove(best_i)
+            print(k, best_degree)
+            self.best_degrees.append(best_degree)
+            if not degree_dic[best_degree]:
+                del degree_dic[best_degree]
+
+        assert len(self.a_ind) == self.n_nodes
+        fig, ax = plt.subplots()
+        ax.plot(self.best_degrees)
+        fig.show()
+        return best_X
+
+    def execute(self, min_iters, min_iters_since_change=1, starting_iters=0, starting_iters_since_change=0, **pars):
+        X, x = self.compute_x()
+        no_a_ind = set(np.arange(self.N)).difference(self.a_ind)
+
+        if self.n_nodes > 1:
+            i = starting_iters
+            iters_since_change = starting_iters_since_change
+            # weights = self.seq / np.sum(self.seq)
+            while i < min_iters or iters_since_change < min_iters_since_change:
+                # new_ind = np.random.randint(self.N)
+                # while new_ind in self.a_ind:
+                #     new_ind = np.random.choice(self.N, replace=False, p=weights)
+
+                ind_untargeted = np.random.choice(list(no_a_ind))
+                ind_targeted = np.random.choice(list(self.a_ind))
+                self.a_ind.remove(ind_targeted)
+                self.a_ind.add(ind_untargeted)
+
+                new_X, new_x = self.compute_x()
+
+                if new_X > X:
+                    X, x = new_X, new_x
+                    no_a_ind.remove(ind_untargeted)
+                    no_a_ind.add(ind_targeted)
+                    iters_since_change = 0
+                else:
+                    self.a_ind.remove(ind_untargeted)
+                    self.a_ind.add(ind_targeted)
+                    iters_since_change += 1
+
+                i += 1
+
+            print("Hill climbing finished after {:.0f} iterations. Iterations since last change: {:.0f}".format(i,
+                                                                                                            iters_since_change))
+        else:
+            print("Discrete setting with a single node, performing greedy algorithm.")
+            max_j = -1
+            for j in range(self.N):
+                self.a_ind = set([j])
+                new_X, new_x = self.compute_x()
+                if new_X > X:
+                    X, x = new_X, new_x
+                    max_j = j
+            self.a_ind = set([max_j])
+            i, iters_since_change = 999999, 999999
+
+        self.produce_a_vec()
+        return self.a[:, 0], X, x, i, iters_since_change
+
+    def compute_x(self, produce_a_vec=True):
+        if produce_a_vec:
+            self.produce_a_vec()
+        inv_Lab = np.linalg.inv(self.L + np.diagflat(self.a + self.b))
+        x = np.array(np.matmul(inv_Lab, self.a).T)
+        X = np.sum(x) / self.N
+
+        return X, x[0, :]
+

networks.py

+import warnings
+import random
+import numpy as np
+import networkx as nx
+import pandas as pd
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+from python.utils import PLOTS_PATH
+from python.gascent import GA, DebugGA
+from python.hillclimb import HillClimber
+
+import os
+
+dirname = os.path.dirname(__file__)
+EMAIL_PATH = os.path.join(dirname, "datasets/ia-email-univ/email.txt")
+ADVOGATO_PATH = os.path.join(dirname, "datasets/soc-advogato/soc-advogato.edges")
+CLASSROOM_PATH = os.path.join(dirname, "datasets/classroom/Friendship-network_data_2013.csv")
+
+mpl.rcParams['text.usetex'] = True
+
+
+def get_highest_nodes(seq, n, end):
+    if end == "high":
+        sorted_indices = np.argsort(seq)[::-1]
+        d_star = seq[sorted_indices][n - 1]
+        chosen_indices = sorted_indices[seq[sorted_indices] > d_star]
+    elif end == "low":
+        sorted_indices = np.argsort(seq)
+        d_star = seq[sorted_indices][n - 1]
+        chosen_indices = sorted_indices[seq[sorted_indices] < d_star]
+    else:
+        raise Exception("end {:s} not understood.".format(end))
+
+    chosen_indices = np.append(chosen_indices, random.sample(
+        list(sorted_indices[seq[sorted_indices] == d_star]), n - len(chosen_indices)))
+
+    return chosen_indices
+
+
+def build_graph(network_class, seed, **pars):
+    if "N" in pars:
+        N = pars['N']
+    else:
+        raise Exception("Number of nodes N not specified!")
+
+    net_pars = {}
+    if network_class == "complete":
+        G = nx.complete_graph(N)
+    elif network_class == "bipartite":
+        G = nx.complete_bipartite_graph(int(N * pars['rho']), int(N * (1 - pars['rho'])))
+    elif network_class == 'ba':
+        G = nx.barabasi_albert_graph(N, pars['m'], seed=seed)
+    elif network_class == 'conf':
+        fails = 0
+        while True:
+            try:
+                # seq = np.array(nx.utils.powerlaw_sequence(N, pars['gamma'], seed=seed), dtype=int)
+
+                # nx.utils.zipf_rv(alpha=pars['gamma'], xmin=pars['dmin'], seed=seed)
+                # print([k for k in np.arange(np.sum(seq == 1))])
+
+                seq = [nx.utils.zipf_rv(alpha=pars['gamma'], xmin=pars['dmin'], seed=seed + k) for k in range(N)]
+                print(max(seq), np.mean(seq), np.unique(seq, return_counts=True))
+                # G = nx.random_degree_sequence_graph(seq)  # , tries=50)
+                G = nx.configuration_model(seq, seed=seed)
+                if not nx.is_connected(G):
+                    raise nx.exception.NetworkXError("Graph not connected")
+                break
+            except (nx.exception.NetworkXError, nx.exception.NetworkXUnfeasible) as err:
+                print(err)
+                seed += 100
+                fails += 1
+        G = nx.Graph(G)  # remove parallel edges
+        G.remove_edges_from(nx.selfloop_edges(G))
+        net_pars["fails"] = fails
+    elif network_class == "regular":
+        G = nx.random_regular_graph(pars['kdegree'], N, seed=seed)
+        # G = nx.fast_gnp_random_graph(pars['N'], 2 * pars['kdegree'] / (pars['N'] - 1))
+    elif network_class == "small-world":
+        G = nx.connected_watts_strogatz_graph(N, pars['kdegree'], pars['p'], seed=seed)
+    elif network_class == "email":
+        edges = pd.read_csv(EMAIL_PATH, sep='\s+').values.tolist()
+        G = nx.Graph(edges)
+        for (u, v) in nx.selfloop_edges(G):
+            G.remove_edge(u, v)
+        # print("N of nodes: " + str(G.number_of_nodes()) + ", N of edges: " + str(G.number_of_edges()))
+    elif network_class == "advogato":
+        G = nx.read_weighted_edgelist(ADVOGATO_PATH, comments="%")  # , create_using=nx.DiGraph)
+        for i, (u, v) in enumerate(nx.selfloop_edges(G)):
+            G.remove_edge(u, v)
+        G = G.subgraph(max(nx.connected_components(G), key=len))
+    elif network_class == "classroom":
+        G = nx.read_edgelist(CLASSROOM_PATH)  # , create_using=nx.DiGraph)
+        G = G.subgraph(max(nx.connected_components(G), key=len))
+    elif network_class == "random":
+        G = nx.fast_gnp_random_graph(N, pars['p'], seed=seed)
+    else:
+        raise Exception("Network class '" + network_class + "' not understood.")
+
+    assert nx.is_connected(G), "Graph not connected"
+    return G, net_pars
+
+
+def compute_gain_nnodes(net, controller="A", **pars):
+    if controller == "A":
+        budget = net.mean_a * net.N
+    elif controller == "B":
+        budget = net.mean_b * net.N
+    else:
+        raise Exception("Controller {:s} not understood.".format(controller))
+
+    if 'N_allocations_a' in pars and controller == "A":
+        if type(pars['N_allocations_a']) == int:
+            n_nodes = pars['N_allocations_a']
+        elif pars['N_allocations_a'] in ("full_shielding", "over_shielding"):
+            n_neighbours = np.squeeze(np.sum(net.W[net.b[:, 0] > 0], axis=0))
+            mask_Nb = np.logical_and(n_neighbours > 0, net.b[:, 0] == 0)
+            n_nodes = np.sum(mask_Nb)
+            if pars['N_allocations_a'] == "over_shielding":
+                n_nodes = int(n_nodes * 1.5)
+        else:
+            raise Exception("Parameter N_allocations_a not understood.")
+        gain = budget / n_nodes
+        print("Allocations from the A-controller:", n_nodes, "gain:", gain)
+    elif 'gain' in pars:
+        n_nodes = int(budget / pars['gain'])
+        if 'gain_adjustment' in pars:
+            if pars['gain_adjustment']:
+                gain = budget / n_nodes
+                print("Gain adjusted to the new value", gain)
+            else:
+                gain = pars['gain']
+        else:
+            warnings.warn("Gain adjustment not defined for these experiments (!)")
+            gain = pars['gain']
+    elif 'N_allocations' in pars:
+        n_nodes = pars['N_allocations']
+        gain = budget / n_nodes
+    else:
+        raise Exception("Discrete allocations not completely defined.")
+
+    return gain, n_nodes
+
+
+# TODO: THERE SHOULD BE AN ABSTRACT NETWORK CLASS AND CHILDREN CORRESPONDING TO SPECIFIC NETWORK TYPES
+class Network:
+    def __init__(self, network_class, seed=5, adjlist=None, node_order=None, **pars):
+        if adjlist is None:
+            self.G, self.construction_graph_data = build_graph(network_class, seed, **pars)
+        else:
+            self.G = nx.parse_adjlist(adjlist)
+        if network_class != "conf":
+            self.W = nx.adjacency_matrix(self.G)  # , range(1, pars['N']+1))
+        else:
+            # print(type(node_order), len(node_order), node_order)
+            # self.W = nx.adjacency_matrix(self.G, node_order)
+            # print("nodes", self.G.nodes)
+            # print("edges", self.G.edges)
+            self.W = nx.adjacency_matrix(self.G)  # , nodelist=np.arange(pars['N']))
+        self.L = np.diagflat(np.sum(self.W, axis=0)) - self.W
+        self.N = self.W.shape[0]
+        if "N" in pars:
+            assert self.N == pars['N'], "{:d} != {:d}".format(self.N, pars['N'])
+        self.seq = np.squeeze(np.array(np.sum(self.W, axis=1)))
+        self.sort_ind = None
+        if network_class != "email":
+            print(network_class, " network with mean degree", np.mean(self.seq), "max degree", np.max(self.seq))
+
+        self.sample = None
+        self.seed = seed
+        self.mean_b = None
+        self.b = None
+
+        self.mean_a = None
+        self.a = None
+        self.ga = None
+        self.x, self.X = None, None
+
+        self.a_disc, self.x_disc = None, None
+        self.hc = None
+
+    def allocate_a(self, a_vec, type_of_experiment):
+        assert type(a_vec) == np.ndarray
+        assert len(a_vec) == self.N
+        assert np.isclose(np.mean(a_vec), self.mean_a), "{:f} {:f}".format(np.mean(a_vec), self.mean_a)
+        # if len(a_vec.shape) == 1:
+        #     a_vec = a_vec[:, None]
+        if type_of_experiment == "GA":
+            self.a = a_vec
+        elif type_of_experiment == "HillClimb":
+            self.a_disc = a_vec
+        else:
+            raise Exception("Type of experiment {:s} not understood.".format(type_of_experiment))
+
+    def allocate_b(self, allocation_type, **pars):  # 'b' is the mean allocation given by the B-controller
+        assert self.mean_b is not None
+        random.seed(self.seed)
+        np.random.seed(self.seed)
+
+        if type(allocation_type) == np.ndarray:
+            self.b = allocation_type
+            assert len(self.b) == self.N
+            assert np.isclose(np.mean(self.b), self.mean_b), "{:f} {:f}".format(np.mean(self.b), self.mean_b)
+            if len(self.b.shape) == 1:
+                self.b = self.b[:, None]
+        else:
+            if allocation_type in ("linear", "powerlaw", "sin", "Amajor", "random", "randn"):
+                self.allocate_b_wrt_func(allocation_type, **pars)
+            elif allocation_type == 'discrete':
+                if 'gain' in pars:
+                    n_nodes = int(self.mean_b * self.N / pars['gain'])
+                elif 'N_allocations' in pars:
+                    n_nodes = pars['N_allocations']
+                else:
+                    raise Exception("Discrete allocations not completely defined.")
+                self.allocate_b_discrete(n_nodes=n_nodes, end=pars['end'])
+            elif allocation_type == "randlin":
+                self.allocate_uniform_b_with_perturbations(pars['n_perturbations'], pars['end'])
+
+            else:
+                raise Exception("B allocation '" + str(allocation_type) + "' not understood.")
+            if 'share_of_uniform' in pars:
+                self.partially_allocate_b_uniformly(pars['share_of_uniform'])
+
+            assert np.isclose(np.mean(self.b), 1), str(np.mean(self.b)) + " != 1"
+            assert np.all(self.b >= 0), str(self.b[self.b < 0])
+            self.b = self.b * self.mean_b
+
+    def allocate_b_wrt_func(self, allocation_type, beta=None, **pars):
+        assert beta is not None or allocation_type in ("random", "randn")
+
+        def beta_allocation(k):
+            if allocation_type == "linear":
+                return k * beta
+            if allocation_type == "powerlaw":
+                return k ** beta
+            elif allocation_type == "sin":
+                assert 0 < beta <= 1
+                return 1 + beta * np.sin(k / np.max(self.seq) * 8 * np.pi)
+            elif allocation_type == "Amajor":
+                assert 0 < beta <= 1
+                return 1 + beta * (np.sin(k / np.max(self.seq) * 8.8 * np.pi) + np.sin(
+                    k / np.max(self.seq) * 11 * np.pi) + np.sin(k / np.max(self.seq) * 13.2 * np.pi)) / 3
+            elif allocation_type == "random":
+                return random.random()
+            elif allocation_type == "randn":
+                allocation = 1 + np.random.rand()
+            elif allocation_type == "randlin":
+                allocation = 1 + 0.005 * np.abs(k - beta) * np.random.randn()
+            else:
+                raise Exception("B allocation '" + str(allocation_type) + "' not understood.")
+            if allocation < 0:
+                allocation = 0
+            elif allocation > 2:
+                allocation = 2
+            return allocation
+
+        b = np.zeros((self.N, 1))
+        random.seed(self.seed)
+        for j, k in enumerate(self.seq):
+            b[j] = beta_allocation(k)
+
+        # Re-scaling
+        if allocation_type in ("linear",):  # "randlin"):
+            # Shift (keeping slope)
+            self.b = b + (self.N - np.sum(b)) / self.N
+
+            while np.any(self.b < 0):
+                broken_nodes = self.b <= 0
+                self.b += (self.N - np.sum(self.b[self.b > 0])) / np.sum(self.b > 0)
+                self.b[broken_nodes] = 0
+        else:
+            # Re-scale
+            self.b = b / np.sum(b) * self.N
+
+    def allocate_uniform_b_with_perturbations(self, n_perturbations, degree):
+        if type(degree) == "string":
+            if degree in ("low, high"):
+                indices = get_highest_nodes(self.seq, 2 * n_perturbations, degree)
+            else:
+                raise ValueError("value not understood: " + degree)
+        elif type(degree) == int:
+            indices = np.arange(self.N)[self.seq == degree]
+        else:
+            ValueError("Error with degree of type: " + str(type(degree)))
+        assert 2 * n_perturbations <= len(indices), "{:d} {:d}".format(n_perturbations, len(indices))
+
+        self.b = np.ones((self.N, 1))
+        for j in range(n_perturbations):
+            i1, i2 = np.random.choice(indices, 2, replace=False)
+            indices = indices[np.logical_not(np.logical_or(indices == i1, indices == i2))]
+
+            perturbation = 1  # 0.3 * np.random.randn()
+            assert perturbation <= 1, "Perturbation too big: " + str(perturbation)
+            self.b[i1] += perturbation
+            self.b[i2] -= perturbation
+
+    def partially_allocate_b_uniformly(self, share):
+        self.sample = np.random.choice(self.N, int(share * self.N), replace=False)
+        bs = np.sum(self.b[self.sample])
+        self.b[self.sample] = bs / int(share * self.N)
+
+    def allocate_b_discrete(self, n_nodes, end):  # end= {'low' degree, 'high' degree, 'middle' = random}
+        b = np.zeros((self.N, 1))
+
+        if end == "middle":
+            b[random.sample(range(self.N), n_nodes)] = 1
+        else:
+            if end == "high":
+                sorted_indices = np.argsort(self.seq)[::-1]
+                d_star = self.seq[sorted_indices][n_nodes - 1]
+                chosen_indices = sorted_indices[self.seq[sorted_indices] > d_star]
+            elif end == "low":
+                sorted_indices = np.argsort(self.seq)
+                d_star = self.seq[sorted_indices][n_nodes - 1]
+                chosen_indices = sorted_indices[self.seq[sorted_indices] < d_star]
+
+            if end in ("high", "low"):
+                chosen_indices = np.append(chosen_indices, random.sample(
+                    list(sorted_indices[self.seq[sorted_indices] == d_star]), n_nodes - len(chosen_indices)))
+            elif end == "periphery":
+                max_neigh_degs = np.zeros(self.N)
+                for i in range(self.N):
+                    max_neigh_degs[i] = np.max(self.W[i] * self.seq)
+                sorted_ind = np.argsort(max_neigh_degs)
+                # sorted_neigh_degs = max_neigh_degs[sorted_ind]
+                chosen_indices = sorted_ind[:n_nodes]
+                # print("degree of nodes", self.seq[chosen_indices])
+                print("max degree of neighbours", max_neigh_degs[chosen_indices])
+                d_star = np.max(self.seq[chosen_indices])
+            else:
+                raise Exception("end {:s} not understood.".format(end))
+
+            b[chosen_indices] = 1
+            print("Discrete b allocation with {:d} nodes, cut degree {:d} ({:s})".format(n_nodes, d_star, end))
+
+        self.b = b / np.sum(b) * self.N
+
+    def gradient_ascent(self, debug=False, **pars):
+        if 'seed' not in pars:
+            pars['seed'] = self.seed
+
+        if not debug:
+            self.ga = GA(self.W, self.b, a_budget=self.mean_a * self.N, **pars)
+        else:
+            if pars['initial_a'] == "random_disc":
+                gain, n_nodes = compute_gain_nnodes(self, **pars)
+                a = np.zeros(self.N)
+                np.random.seed(self.seed)
+                a[np.random.choice(self.N, n_nodes, replace=False)] = gain
+                pars['initial_a'] = a
+            self.ga = DebugGA(self.W, self.b, a_budget=self.mean_a * self.N, **pars)
+        self.a, self.X, self.x, i = self.ga.execute()
+        return self.a, self.X, self.ga.conv, i
+
+    def hill_climbing(self, n_nodes, gain, min_iters, **pars):
+        self.hc = HillClimber(self.W, self.b, n_nodes, gain, **pars)
+        self.a_disc, X, self.x_disc, iters, i_since_change = self.hc.execute(min_iters, **pars)
+        assert np.isclose(np.mean(self.a_disc), self.mean_a)
+        return self.a_disc, X, self.x_disc, iters, i_since_change
+
+    def sort_by_degree(self, type_of_experiment="GA"):
+        if self.sort_ind is None:
+            self.sort_ind = np.argsort(self.seq)
+        if type_of_experiment == "GA":
+            return self.seq[self.sort_ind], self.a[self.sort_ind], self.b[self.sort_ind, 0], self.x[self.sort_ind]
+        elif type_of_experiment == "HillClimb":
+            return self.seq[self.sort_ind], self.a_disc[self.sort_ind], self.b[self.sort_ind, 0], self.x_disc[
+                self.sort_ind]
+        else:
+            return self.seq[self.sort_ind], self.sort_ind
+
+    def get_edge_list(self):
+        return list(nx.generate_adjlist(self.G))
+
+    @staticmethod
+    def extract_network_pars(network_class, **pars):
+        net_pars = {}
+        if network_class == 'ba':
+            net_pars["m"] = pars['m']
+        elif network_class == 'conf':
+            net_pars["gamma"] = pars['gamma']
+        elif network_class == "regular":
+            net_pars["kdegree"] = pars["kdegree"]
+        elif network_class == "small-world":
+            net_pars["kdegree"] = pars['kdegree']
+            net_pars["p"] = pars['p']
+        elif network_class == "random":
+            net_pars["p"] = pars['p']
+        elif network_class != "email":
+            raise Exception("Network class '" + network_class + "' not understood.")
+
+        return net_pars
+
+
+def create_regular_network(pars):
+    G = nx.random_regular_graph(pars['kdegree'], pars['N'])
+    # G = nx.fast_gnp_random_graph(pars['N'], 2 * pars['kdegree'] / (pars['N'] - 1))
+    if not nx.is_connected(G):
+        raise Exception("Graph not connected")
+    W = nx.adjacency_matrix(G)
+    L = np.diagflat(np.sum(W, axis=0)) - W
+    N = W.shape[0]
+
+    return W, L, N
+
+
+if __name__ == "__main__":
+    pars = {
+        'N': 1133,
+        'p': 0.01,  # For random graph
+        'gamma': 3, 'dmin': 5,  # For scale-free conf model
+        'm': 5,  # For Barabasi-Albert network
+    }
+
+    # network_type = "email"
+    # network_type = "advogato"
+    network_type = "classroom"
+
+    # network_type = "random"
+    # network_type = "bamean"
+    # network_type = "conf"
+
+    if network_type in ("random", "ba", "conf"):
+        n_seeds = 15
+    else:
+        n_seeds = 1
+
+    if network_type == "advogato":
+        pars['N'] = 5054
+    elif network_type == "classroom":
+        pars["N"] = 128
+
+    for k in range(n_seeds):
+        seed = 6 + k
+        net = Network(network_type, seed=seed, **pars)
+        print(net.N)
+        print("degrees (max, mean, seq)", max(net.seq), np.mean(net.seq), np.unique(net.seq, return_counts=True))
+        print("assortativity", nx.degree_assortativity_coefficient(net.G))
+
+    if True:
+        if False:
+            fig, ax = plt.subplots(figsize=(1.618 * 2.5, 2.5))
+            ax.hist(net.seq, bins=range(int(np.max(net.seq)) + 1), color="purple")
+            ax.set_yscale("log")
+            ax.set(xlabel="Node degree", ylabel="Count")
+            ax.margins(x=0.0)
+
+            fig.tight_layout()
+            fig.savefig(PLOTS_PATH + network_type + "_degree_dist.eps")
+        else:
+            print("Directed?", nx.is_directed(net.G))
+            fig, ax = plt.subplots(figsize=(3, 2.5))
+            ax.bar(range(net.N), sorted(net.seq)[::-1], color="purple")
+            ax.set_yscale("log")
+            ax.set(xlabel=r"$i$ (sorted by $d_{i}$, descending)", ylabel=r"Weighted degree, $d_{i}$")
+            ax.set(xlim=[-10, net.N + 9])
+            fig.tight_layout()
+            fig.savefig(PLOTS_PATH + network_type + "_degree_dist_nodes.png")
+
+        fig.show()
+    else:
+        fig, ax = plt.subplots(figsize=(10, 5))
+
+        net.mean_b = 1
+        net.allocate_b("randlin", n_perturbations=10, end=1)
+        # print(np.unique())
+
+        ax.scatter(net.seq, net.b)
+
+        ax.set(ylim=0)
+        fig.show()

requirements.txt

+cycler==0.10.0
+decorator==4.4.2
+kiwisolver==1.2.0
+matplotlib==3.2.2
+mpmath==1.1.0
+networkx==2.4
+numpy==1.19.0
+pandas==1.0.5
+pyparsing==2.4.7
+python-dateutil==2.8.1
+pytz==2020.1
+scipy==1.5.1
+seaborn==0.10.1
+six==1.15.0
+sympy==1.6.2

utils.py

+import numpy as np
+import networkx as nx
+
+import os
+dirname = os.path.dirname(__file__)
+DATA_PATH = os.path.join(dirname, 'simulation_data/')
+PLOTS_PATH = os.path.join(dirname, 'plots/')
+
+
+def std_errors(array, axis=0, confidence=False, nan=False):
+    if type(axis) in (tuple, list):
+        n_dim_averaged = sum(array.shape[ax] for ax in axis)
+    else:
+        n_dim_averaged = array.shape[axis]
+
+    if nan:
+        errors = np.nanstd(array, axis=axis) / np.sqrt(n_dim_averaged)
+    else:
+        errors = np.std(array, axis=axis) / np.sqrt(n_dim_averaged)
+
+    if confidence:
+        return 2 * errors
+    else:
+        return errors
+
+
+def create_range(pars, level='axis'):
+    if level == "axis":
+        if pars[level + "par"] in ("rho", "beta") and level + "array" not in pars:
+            for el in np.linspace(pars[level + "min"], pars[level + "max"], pars[level + 'length']):
+                pars[pars[level + 'par']] = el
+                yield el
+        elif pars[level + 'par'] == "k1":
+            for rho in np.linspace(1 / pars['N'], (pars['N'] - 1) / pars['N'], pars['N'] - 1):
+                num = (pars['kdegree'] - rho * pars['k0']) / (1 - rho)
+                if num < 28 + 7 * pars['k0'] and np.isclose(np.mod(num, 1), 0):
+                    seq = [pars['k0'] for _ in range(int(pars['N'] * rho))] + [int(num) for _ in
+                                                                               range(int(pars['N'] * (1 - rho)))]
+                    if nx.is_graphical(seq):
+                        pars['k1'] = int(num)
+                        pars['rho'] = rho
+                        if 'betamul' in list(pars):
+                            pars['alpha'] = rho + (1 - rho) * (1 - pars['betamul'])
+                            pars['beta'] = pars['betamul'] * rho
+                        yield (int(num), rho)
+        elif level + "array" in pars:
+            for el in pars[level + 'array']:
+                pars[pars[level + 'par']] = el
+                yield el
+        else:
+            raise Exception("Axis parameter not understood: " + pars[level + "par"])
+    elif level == "lines":
+        if level + "array" in pars:
+            for el in pars[level + 'array']:
+                if pars[level + 'par'] is not None:
+                    pars[pars[level + 'par']] = el
+                yield el
+        else:
+            raise Exception("Lines parameter '{:s}' not understood.".format(pars[level + 'par']))
+    elif level == "axes" or level == "multiaxes":
+        if pars[level + "par"] == "ab":
+            for el in np.power(10.0, np.array([-1, 0, 1])):
+                pars['a'], pars['b'] = el, el
+                yield el
+        elif level + "array" in pars:
+            for el in pars[level + 'array']:
+                if pars[level + 'par'] is not None:
+                    pars[pars[level + 'par']] = el
+                yield el
+        else:
+            raise Exception("Axes par '{:s}' not understood".format(pars[level + "par"]))
+    else:
+        raise Exception("Level '{:s}' not understood.".format(level))
+