Gohモデル。

Gohモデル、任意の次数分布をもつスケールフリーネットワークを作る方法。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>





int main(){


	int N=1000,i;
	int k_mean=3;
	double* omega;
	double alpha=2.3;
	double* p;
	double* q;
	int* a;
	
	int* h;
	int* l;
	int j;
	double* pcum;
	int* x;
	double sum=0;
	int n;
	int lj;
	int hj;
	int counta;
	int indexi,indexj;
	int u2;
	double u1;
	
	int *out;
	int *in;

	int L=k_mean*N/2;

	omega=(double*)calloc(5*N,sizeof(double*));

	p=(double*)calloc(5*N,sizeof(double*));

	q=(double*)calloc(5*N,sizeof(double*));
	a=(int*)calloc(5*N,sizeof(int*));
	x=(int*)calloc(5*N,sizeof(int*));
	h=(int*)calloc(5*N,sizeof(int*));
	l=(int*)calloc(5*N,sizeof(int*));

	out=(int*)calloc(2*L,sizeof(int*));

	in=(int*)calloc(2*L,sizeof(int*));

	for(i=0;i<=N-1;i++){

		omega[i]=pow(i+1.0,-alpha);			
		x[i]=i+1;
		sum=sum+omega[i];
	}
	
	for(i=0;i<=N-1;i++){

		p[i]=omega[i]/sum;
		q[i]=(double)N*p[i];		
		a[i]=i+1;
	}
	lj=0;
	hj=0;

	for(i=0;i<=N-1;i++){
		if(q[i]>1){	
			h[hj]=i;
			hj++;	

			//printf("h%d\n",i);
		}else{	
			l[lj]=i;
			lj++;	
			//printf("l%d\n",i);
		}
	}	

	indexi=0;
	indexj=0;

	for(i=0;i<=hj-1;i++){
		//printf("%d %d\n",h[i],hj);
	}	

	


	while(indexi<hj){
		i=h[indexi];
		j=l[indexj];

		a[j]=x[i];
		q[i]=q[i]-(1-q[j]);

		if(q[i]<1){
			indexi++;
			if(l[lj]!=i){
				l[lj]=i;
				lj++;
			}
		}


		indexj++;	
		

	}


	
	srand(1);

	int node;
	int node2;
	int i12=0;
	int T=1;

	while(i12<=L){
		u1=(double)rand()/RAND_MAX;	
		u2=rand()%(N);	
		if(q[u2] >= u1){
			node=x[u2];

		}else{
			node=a[u2];
		}	
		

		u1=(double)rand()/RAND_MAX;	
		u2=rand()%(N);	
		if(q[u2] >= u1){
			node2=x[u2];

		}else{
			node2=a[u2];
		}	
	
		T=1;
		for(i=0;i<=i12-1;i++){
			if(out[i]==node && in[i]==node2){

				T=0;
			}
				

		}	

		if(node!=node2 && T){
			out[i12]=node;
			in[i12]=node2;
			i12++;
		}

		
	}	
	for(i=0;i<=L-1;i++){
		printf("%d,%d\n",out[i],in[i]);
	}

}

下のコードのが速い。
c++を利用

#include<stdio.h>
#include<math.h>
#include<string>
#include<string.h>
#include<vector>
#include <iostream>
#include<set>
#include<map>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<stdlib.h>

int main(){
        using namespace std;
        const gsl_rng_type *T;
        gsl_rng *r;
        T=gsl_rng_default;
        r=gsl_rng_alloc(T);
        gsl_rng_set(r,1);

        int *out1;
        int *in1;
        double *omega;
        double *p;
        double gamma=4.0;
        double alpha=1.0/(gamma-1);
        //double alpha=1.3;
        int i;
        double sumtmp;
        int n123=3000;

        std::map< int, std::set<int> > adjency_list;
        double meank=(gamma-1)/((gamma-1)-1.0);

        out1=(int*)calloc((int)(10*n123*meank+1),sizeof(int));
        in1=(int*)calloc((int)(10*n123*meank+1),sizeof(int));

        omega=(double*)calloc(n123,sizeof(double));
        p=(double*)calloc(n123,sizeof(double));

        sumtmp=0;
        for(i=0;i<n123;i++){
                omega[i]=pow(i+1,-alpha);
                sumtmp=sumtmp+omega[i];

        }


        for(i=0;i<n123;i++){

                p[i]=omega[i]/sumtmp;
                //cout <<omega[i]<<" "<<p[i] << endl;
        }


        gsl_ran_discrete_t* t=gsl_ran_discrete_preproc(n123,p);
        int i1;
        int i2;
        i=0;
        while(i<meank*n123/2){

                i1= gsl_ran_discrete(r,t) ;
                i2= gsl_ran_discrete(r,t) ;

                //if(adjency_list[i1].find(i2)==adjency_list[i1].end() && adje
//y_list[i2].find(i1)==adjency_list[i2].end() && i1!=i2){

                if(adjency_list[i1].find(i2)==adjency_list[i1].end() && i1!=i2){

                        out1[i]=i1;
                        in1[i]=i2;

                        adjency_list[i1].insert(i2);


                        i++;
                }

        }

        for(i=0;i<meank*n123/2;i++){
                cout << out1[i] <<" "<<in1[i] << endl;

        }


}