//Create by Wubo
//Revise on 08/25/06 03:32

#include "d:\scilab-4.0\routines\stack-c.h"
void sis __PARAMS((int *a, int n, int m, int *b, int *nodenum,int *step,int *n0, double *PB,double *inf,int *infnode));

int intsis(fname) 
char *fname;
{
    static int l1, m1, n1, l2, m2, n2, l3,m3,n3,l4,l5,l6,l7,l8;
    static int minlhs=0, maxlhs=2, minrhs=6, maxrhs=6;   //why the minlhs must equal zero??

    /* Check number of inputs (rhs=2) and outputs (lhs=1) */
    CheckRhs(minrhs,maxrhs) ;
    CheckLhs(minlhs,maxlhs) ;
    //m2=10;
    GetRhsVar(1, "i", &m1, &n1, &l1);   //the following sentence must be wrote together!!!! 
    GetRhsVar(2, "i", &m1, &n1, &l2);  //*b 
    GetRhsVar(3, "i", &m1, &m1, &l3);  //*nodenum
    GetRhsVar(4, "i", &m1, &m1, &l4);  //*step
    GetRhsVar(5, "i", &m1, &m1, &l5);  //*n0;
    GetRhsVar(6, "d", &m1, &m1, &l6);  //*PB
    m2=*istk(l4);
    m3=*istk(l3);
    CreateVar(7, "d", &m1, &m2, &l7); 
    CreateVar(8, "i", &m1, &m3, &l8); 
    sis(istk(l1), m1, n1, istk(l2),istk(l3),istk(l4),istk(l5),stk(l6),stk(l7),istk(l8));
    /*  Return C (3)  */
    LhsVar(1) = 7;
    LhsVar(2) = 8;
    return 0;
}


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

int N;
int T;
int N0;
double P;


int      **node;
int      *state;        //state of the node {S==0,I==1,R==2}
int      *temp;
int      *nodetime;     //account the time of each node
double   *inf;          //fraction of I 
double   *avedeg;       //ave deg of infected nodes    
int      *deg; 
int *inftimes;
int *infname;




void input(int *head,int *tail,int edge,int N)
{
    int aa,bb;
    int i,j;
    int e=0;
    for (i=0;i<N;i++) {
        deg[i]=0;
        for (j=0;j<2000;j++) {
            node[i][j]=0;
        }
    }
    while (e<edge) {
        aa=head[e]-1;
        bb=tail[e]-1;
        node[aa][deg[aa]]=bb;
        node[bb][deg[bb]]=aa;
        deg[aa]++;
        deg[bb]++;
        e++;
    }
}


void initial(int T,int N)
{
    int i,temp2,num;
    for (i=0;i<T;i++) {
        inf[i]=0;
        avedeg[i]=0;
    }
    for (i=0;i<N;i++) {
        state[i]=temp[i]=0;
    }
    num=0;
    while (num<N0) {
        temp2=rand()%N;
        if (state[temp2]==0) {
            state[temp2]=temp[temp2]=1;
            num++;
        }
    }

}

void infection(int t,int N)
{
    int i,j;
    int ti=0;
    int index,tempnod;
    double frand;
    for (i=0;i<N;i++) {
        state[i]=temp[i];
    }
    for (i=0;i<N;i++) {
        temp[i]=0;
    }

    for (i=0;i<N;i++) {
        if (state[i]==1) {
            for (j=1;j<=node[i][0];j++) {

                frand=(double)rand()/RAND_MAX;
                if (frand<P) {
                    temp[j]=1;
                    inftimes[j]++;
                    ti++;
                    avedeg[t]+=deg[j];
                }

            }
        }
    }

    index=0;

    for (i=0;i<N;i++)
        infname[i]=-1;


    for (i=0;i<N;i++)
        if (state[i]==1) {
            infname[index]=i;
//			 sciprint("%i \n",i);
            index++;
            tempnod=i;
        }
    for (i=index;i<N;i++)
        infname[i]=tempnod;



    avedeg[t]/=ti;
    for (i=0;i<N;i++) {
        state[i]=0;
    }
}

void  tongji(int av,int t,int N)
{
    int i;
    for (i=0;i<N;i++) {
        if (temp[i]==1) {
            inf[t]++;
        }
    }
    inf[t]/=N;

}








void sis(a,n,m,b,nodenum,step,n0,PB,infm,infnode)
int a[],b[],nodenum[],step[],n0[];
double infm[],PB[];
int n,m;
int infnode[];
{
    int i,j;

    N=*(nodenum);
    T=*(step);
    N0=*(n0);
    P=*(PB);

    node=malloc(N*sizeof(int *));
    for (i=0;i<N;i++)
        node[i]=malloc(2000*sizeof(int));   //generate node[N][2000]

    state=(int *) malloc(N * sizeof(int));
    temp=(int *)malloc(N*sizeof(int));
    inf=(double *)malloc(T*sizeof(double));
    avedeg=(double *)malloc(T*sizeof(double));
    inftimes=(int *)malloc(N*sizeof(int));
    deg=(int *)malloc(N*sizeof(int));
    infname=(int *)malloc(N*sizeof(int));

    srand((unsigned)time(NULL));
    input(a,b,m,N);

    for (i=0;i<N;i++) {
        inftimes[i]=0;
    }


    for (j=0;j<1;j++) {
        initial(T,N);



        for (i=0;i<T;i++) {
            infection(i,N);
            tongji(j,i,N);
        }

        infm[0]=0;
        for (i=1;i<T;i++) {
            infm[i]=inf[i];
        }
    }



    for (i=0;i<N;i++)
        infnode[i]=infname[i];





    free(node);
    free(state);
    free(temp);
    free(nodetime);
    free(inf);
    free(avedeg);
    free(deg);
    free(inftimes);
    free(infname);


}