📄 simple.cpp
字号:
#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
//#define myabs(x) (((x)>=0)?(x):(-x))
#define ABS(x) (x>0 ? x : -x)
#define min(x,y) (x<y ? x : y)
#define max(x,y) (x>y ? x : y)
#define MAX_WIDTH 10000
/*
* the number of 1s in the integer represented by binary version
* 求解二进制数表示中1的个数
*/
int myabs(int x)
{
if(x>=0)
return x;
else
return -x;
}
unsigned num_1(int b)
{
int b_tmp,num_b,remainder_b;
b_tmp=b;
num_b=0;
while(b_tmp>0)
{
remainder_b=b_tmp%2;
if(remainder_b==1)
num_b+=1;
b_tmp/=2;
}
return num_b;
}
/*
* the probobility of a input vector represented by a integer b
* 用整数b表示的输入向量的概率（二进制数表示中1的个数与输入位数的比）
*/
double pro_b(int b,int n)
{
double b_pro;
b_pro=(double)num_1(b)/((double)n);
return b_pro;
}
/*
* the spatial correction of a input vector represented by a integer b
* 用整数b表示的输入向量的空间相关性
*/
double spa_b(int b,int n)
{
double b_spa;
b_spa=(double)num_1(b)*(n-(double)num_1(b))/ceil((double)(n/2))/floor((double)(n/2));
return b_spa;
}
/*justify whether s is satisfy the restriction */
bool s_satisfy(double p,double s,int n)
{
double s_bound;
/* 大于n*p的最小整数和小于n*p的最大整数 */
/* the minimal integer that bigger than n*p and the max integer that smaller than n*p */
int ceil_np,floor_np;
ceil_np=(int)ceil(n*p);
floor_np=(int)floor(n*p);
if(ceil_np==floor_np)
{
floor_np=(int)(n*p);
ceil_np=floor_np+1;
}
s_bound=ceil_np*(n-ceil_np)*(n*p-floor_np)+floor_np*(n-floor_np)*(ceil_np-n*p);
s_bound=s_bound/(ceil(n/2)*floor(n/2));
if(s_bound>=s && s>=0 && s<=1 && p>=0 && p<=1)
return 1;
else
return 0;
}
/*
* caculate the transition density from cluster k to clastuer j when there are x bits converting
*/
double trans_kj(int n,int x,int k,int j,double lamda)
{
double pro_x,nom_value=0;
int i;
int bottom,top;
bottom=myabs(k-j);
top=min(k+j,2*n-k-j);
for(i=bottom;i<=top;i+=2)
{
nom_value+=pow(lamda,i-bottom)*pow(1-lamda,top-i);
}
pro_x=pow(lamda,x-bottom)*pow(1-lamda,top-x)/nom_value;
return pro_x;
}
/*
* generate a vector j from a initial vector k swiching x bits
*/
bool* vec_gen(int k,int j,int n,int x,bool* initial)
{
/* bit swich from 0 to 1 and that from 1 to 0 */
int bit_10,bit_01;
/* the beginning bit in the vector */
int begin_10,begin_01;
/* output vector */
bool* out_vec;
int i;
out_vec= (bool *) calloc(n+1,sizeof(bool));
for(i=0;i<=n;i++)
{
out_vec[i]=initial[i];
}
bit_01=(x-k+j)/2;
bit_10=(x+k-j)/2;
begin_01=rand()%n;
i=0;
while (i<bit_01)
{
if(initial[begin_01]==0)
{
out_vec[begin_01]=1;
i++;
}
if(begin_01==n-1)
{
begin_01=0;
}
else
{
begin_01++;
}
}
begin_10=rand()%n;
i=0;
while (i<bit_10)
{
if(initial[begin_10]==1)
{
out_vec[begin_10]=0;
i++;
}
if(begin_10==n-1)
{
begin_10=0;
}
else
{
begin_10++;
}
}
return out_vec;
}
/*
* cacluate average transtion density d
*/
double aver_trans(int n,double *trans,double *t,double lamda)
{
double aver_d=0,d_lamda;
int k,j,i;
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
{
d_lamda=0;
for(i=myabs(k-j);i<=min(k+j,2*n-j-k);i+=2)
{
d_lamda+=i*trans_kj(n,i,k,j,lamda);
}
d_lamda/=n;
aver_d+=d_lamda*trans[k*(n+1)+j]*t[k];
}
}
return aver_d;
}
/*
* write the vector to the file
*/
void write_to_file(int n,int i,bool *vec_seq,FILE *infile)
{
int k;
for(k=0;k<n;k++)
{
fprintf(infile,"%d",vec_seq[k]);
}
fprintf(infile,"\n");
}
void seq_gen(int seed,int n,int l,double p0,double d0,double s0)
{
/* the precision of the evluation of t_min */
double precision=1e-6;
/* the bounds of the t_min,*/
double top=0,bottom=1,t_min=0,p_tmp=0,s_tmp=0,d_tmp=0;
int i,j,k;
double sum_pk=0,sum_sk=0;
top=(double)1/((double)(n+1));
bottom=0;
/*evaluate the minimal tk 求最小的tk值*/
printf("calculating t_min\n,precision=%f,top=%f\n",precision,top);
for(k=0;k<=n;k++)
{
sum_pk+=(double)(k)/((double)(n+1));
sum_sk+=spa_b(k,n)/(n+1);
}
while (top-bottom>=precision)// && p_tmp>=0 && s_tmp>=0 && p_tmp<=1 && s_tmp<=1)
{
t_min=(top+bottom)/2;
/* eqution 32,33*/
p_tmp=(p0-(n+1)*t_min*((double)sum_pk)/n)/(1-(n+1)*t_min);
s_tmp=(s0-(n+1)*t_min*sum_sk)/(1-(n+1)*t_min);
/* refresh the bounds of t_min */
if (s_satisfy(p_tmp,s_tmp,n)==1)
{
bottom=t_min;
}
else
{
top=t_min;
}
}
printf("t_min=%f\n",t_min);
/*
* evaluation of t
* 求t的值
*/
/* the vector t,t1;*/
printf("caculating t\n");
double *t;
t= (double *) calloc(n+1,sizeof(double));
for(k=0;k<=n;k++)
{
t[k]=0.0;
}
/* integer bound of n*p0 */
int n1,n2;
/* there's two cases:stmp=0;stmp>0 */
if (0==s_tmp)
{
for(k=1;k<n;k++)
{
t[k]=t_min;
}
t[0]=t_min+(1-(n+1)*t_min)*(1-p0);
t[n]=t_min+(1-(n+1)*t_min)*p0;
}
else
{
if ((int)ceil(n*p0)==(int)floor(n*p0))
{
n2=(int)floor(n*p0);
n1=n2+1;
}
else
{
n1=(int)ceil(n*p0);
n2=(int)floor(n*p0);
}
/* pay attention that n1=n2+1! */
for(k=0;k<=n;k++)
{
t[k]=t_min;
}
t[n1]=t_min+(1-(n+1)*t_min)*(n*p0-n2);
t[n2]=t_min+(1-(n+1)*t_min)*(n1-n*p0);
}
printf(" t caculated\n");
/*
* caculate the transition matrix trans
*/
/* transition matrix and the independent transition matrix,biggest trans */
double *trans;
double *ind_trans;
double *big_trans;
FILE *file_t;
file_t=fopen("file_t","w");
double d_max,d_min,big_dmax;
/* some variables about lamda */
double lamda;
trans=(double *) calloc((int)pow(n+1,2),sizeof(double));
ind_trans=(double *) calloc((int)pow(n+1,2),sizeof(double));
big_trans=(double *) calloc((int)pow(n+1,2),sizeof(double));
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
{
ind_trans[k*(n+1)+j]=t[j];
}
}
/* caculate d_max and d_min*/
d_max=aver_trans(n,ind_trans,t,1);
d_min=aver_trans(n,ind_trans,t,0);
/* caculate the transition density matrix trans*/
/* coefficient to adjust ind_trans */
double co1,co2;
/* the the row sum of big_trans */
double *sum;
double delta=0;
int temp=0;
sum=(double *) calloc(n+1,sizeof(double));
if(d0>d_min && d0<d_max)
{
printf("d0 is middle\n");
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
trans[k*(n+1)+j]=ind_trans[k*(n+1)+j];
}
}
else if(d0<d_min)
{
printf("d0 is small\n");
co1=(d0/d_min+d0/d_max)/2;
for(k=0;k<=n;k++)
for(j=0;j<=n;j++)
{
if(k==j)
trans[k*(n+1)+j]=(1-co1)+co1*ind_trans[k*(n+1)+j];
else
trans[k*(n+1)+j]=co1*ind_trans[k*(n+1)+j];
}
}
else if(d0>d_max)
{
/*
* this arithmetic satisfy (23),but not (22)
*/
printf("d0 is large\n");
for(k=0;k<=n;k++)
{
if(t[temp]<=t[k])
temp=k;
}
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
{
if (j==n-k)
{
big_trans[k*(n+1)+j]=t[temp];
big_trans[k*(n+1)+temp]=t[j];
}
else if(j!=temp)
{
big_trans[k*(n+1)+j]=t[i];
}
sum[j]+=big_trans[k*(n+1)+j]*t[k];
}
}
for(j=0;j<=n;j++)
{
for(k=0;k<=n;k--)
{
big_trans[k*(n+1)+j]=big_trans[k*(n+1)+j]*t[j]*sum[j];
}
}
/*caculate the big_dmax,*/
big_dmax=aver_trans(n,big_trans,t,1);
printf("big_dmax=%d\n",big_dmax);
if(d0>big_dmax)
{
printf("can't generate transition matrix\n");
exit(0);
}
else
{
co2=(1+(d0/d_max)/(big_dmax/d_max))/2;
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
trans[k*(n+1)+j]=co2*big_trans[k*(n+1)+j]+(1-co2)*big_trans[k*(n+1)+j];
}
printf("transition matrix caculated\n");
}
}
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
{
fprintf(file_t,"%f ",trans[k*(n+1)+j]);
}
fprintf(file_t,"\n");
}
free(ind_trans);
free(big_trans);
/* caculate lamda */
top=1;
bottom=0;
precision=0.000001;
while (top-bottom>=precision)
{
lamda=(top+bottom)/2;
d_tmp=aver_trans(n,trans,t,lamda);
if(d_tmp>=1 || d_tmp<=0)
{
printf("can't generate lamda\n");
exit(0);
}
else if(d_tmp>d0)
top=lamda;
else if(d_tmp<d0)
bottom=lamda;
else
{
top=lamda;
bottom=lamda;
}
}
printf("lamda caculated,lamda=%f\n",lamda);
/*
* generate the vector sequence
*/
/*pointing to file that store the generated sequences*/
FILE *infile,*output;
/* seed used to start the generation step and current states identify by "1"the sequence has */
int seq_seed=seed,cur_state=0;
/* 所需l个矢量中各个状态的矢量个数；某状态已经向前trans_num个状态跳转过 */
int *state_num,*trans_num;
/* the number of every transition storing in a matrix */
int *int_trans,*swich_bit;
/*the generated vector */
bool *vec_seq;
bool flag;
/* the swich bits from state k to state j */
int x=0;
vec_seq=(bool *) calloc(n,sizeof(bool));
trans_num=(int *) calloc(n+1,sizeof(int));
state_num=(int *) calloc(n+1,sizeof(int));
int_trans=(int *) calloc((n+1)*(n+1),sizeof(int));
int *x_able_gen;
x_able_gen=(int *) calloc((int)pow(n+1,2),sizeof(int));
swich_bit=(int *) calloc((int)pow(n+1,2),sizeof(int));
/* initializition */
output=fopen("output","w");
for(k=n;k>=0;k--)
{
for(j=n;j>=0;j--)
{
int_trans[k*(n+1)+j]=(int)(ceil(t[k]*trans[k*(n+1)+j]*l));
swich_bit[k*(n+1)+j]=myabs(k-j);
x_able_gen[k*(n+1)+j]=(int)ceil(trans_kj(n,myabs(k-j),k,j,lamda)*int_trans[k*(n+1)+1]);
fprintf(output,"k=%d,j=%d,x_able_gen=%d,int_trans=%ld,trans=%f,swich_bit=%d\n",
k,j,x_able_gen[k*(n+1)+j],int_trans[k*(n+1)+j],trans[k*(n+1)+j],swich_bit[k*(n+1)+j]);
fprintf(file_t,"%d ",int_trans[k*(n+1)+j]);
}
fprintf(file_t,"\n");
state_num[k]=(int)ceil(t[k]*l);
if(rand()%(n+1)==0)
trans_num[k]=n1;
else
trans_num[k]=n2;
printf("%d ",trans_num[k]);
}
for(k=0;k<n;k++)
{
vec_seq[k]=(bool)(seq_seed%2);
seq_seed/=2;
if(1==vec_seq[k])
cur_state++;
}
fprintf(file_t,"initializition finished\n");
i=1;
flag=0;
infile=fopen("sequence","w");
write_to_file(n,0,vec_seq,infile);
while (i<l)
{
j=trans_num[cur_state];
while(int_trans[cur_state*(n+1)+j]<=0 && flag==0)
{
if(j==n)
j=0;
else
j++;
if(j==trans_num[cur_state])
flag=1;
}
if(flag==0 && int_trans[cur_state*(n+1)+j]>0)
{
while(x_able_gen[cur_state*(n+1)+j]<=(int)ceil(l*trans[cur_state*(n+1)+j])-int_trans[cur_state*(n+1)+j]
&& swich_bit[cur_state*(n+1)+j]<min(cur_state+j,2*n-cur_state-j))
{
swich_bit[cur_state*(n+1)+j]+=2;
x_able_gen[cur_state*(n+1)+j]+=(int)ceil(l*t[cur_state]*trans[cur_state*(n+1)+j]*
trans_kj(n,swich_bit[cur_state*(n+1)+j],cur_state,j,lamda));
}
vec_seq=vec_gen(cur_state,j,n,swich_bit[cur_state*(n+1)+j],vec_seq);
write_to_file(n,i,vec_seq,infile);
/* refresh the parametor */
int_trans[cur_state*(n+1)+j]--;
state_num[j]--;
if(trans_num[cur_state]==n)
{
trans_num[cur_state]=0;
}
else
trans_num[cur_state]++;
cur_state=j;
i++;
}
else if(flag==1)
{
printf("can't generate sequence\n");
exit(0);
}
}
for(k=0;k<=n;k++)
{
for(j=0;j<=n;j++)
{
fprintf(file_t,"%d ",int_trans[k*(n+1)+j]);
}
fprintf(file_t,"\n");
}
free(int_trans);
free(swich_bit);
free(state_num);
free(trans_num);
free(vec_seq);
free(x_able_gen);
fclose(infile);
free(infile);
}
int main(int argc, char *argv[])
{
double p0; // average p
double s0; // average s
double d0; //average d
int n; // width
int l; // length
long seed;
FILE * infile;
float ptmp,dtmp,stmp;
if(argc!=3)
{
printf("USAGE: seq_gen inputfile seed\n");exit(0);
}
seed=atoi(argv[2]);
infile=fopen(argv[1],"r");
if(infile==NULL)
{
printf("file can not be opened\n");
exit(0);
}
// process the input parameters
fscanf(infile,"%ld",&n);
fscanf(infile,"%ld",&l);
printf("n=%ld,l=%ld,MAX_WIDTH=%d\n",n,l,MAX_WIDTH);
if(n>MAX_WIDTH)
{
printf("MAX width exceeded\n");exit(0);
}
while(fscanf(infile,"%f %f %f",&ptmp,&dtmp,&stmp)==3)
{
p0=ptmp;
d0=dtmp;
s0=stmp;
}
printf("p0=%f,d0=%f,s0=%d\n",p0,d0,s0);
printf("begin to generate vector sequence!\n");
seq_gen(seed,n,l,p0,d0,s0);
printf("generation finished!\n");
fclose(infile);
return 1;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -