/*
  PROGRAMA PER UN MODEL D'APRENENTATGE HISTORIC-GEOGRAFIC
*/
 
#include "time.h"
#include "stdio.h"
 
#define PC 1       /* 0=PC Convertible   1=PC amb disc dur  2=IBM 3090 */
#define mx0 1      /* dimensio inicial */
#define U0  0      /* valor inicial */
#define mm 4       /* cota superior per la dimensio maxima */
#define NP 100     /* poblacio */
#define MM 15      /* cota superior per l'estat maxim = (1<<mm)-1 */
/*  #define NULL 0  */
 
short int mx=mx0,m[NP],MXz,Mz[NP],NM;
int natal[NP],tanatos=20000,progres=10000,ka,kr,kn;
long int T=0,Ti,F[MM+1][NP],B[NP];
char imprmx=1,neonat[NP],peek();
float P[MM+1][NP],PI[MM+1],st[MM+1],em[MM+1],EMz[NP],a[NP],bz[NP];
FILE *pf; char fnom[35];
float sts[MM+1][NP],ferocitat[MM+1],PR[NP],SR[NP];
float PR0[mm+1],SR0[mm+1];
float IMPA[NP][NP],imp[MM+1][NP/2+1],aT[NP],Tr[NP],K1=100,K2=0,k=1,ki;
 
main()
{randomize();
 presentacio();
 inicialitzar();
 interrupcions();
 valors();
 estadistica();
 impacte();
 /* mostrar(); */
 /* pre_impr(); */
 iteracio();
}
 
randomize()
{void srand(),exit(); time_t time();
 srand((unsigned int)time(NULL)/2);
 #if PC==1
  sprintf(fnom,"C:\\TUROUT\\MODAPRHG\\%lx.OUT",time(NULL));
 #elif PC==0
  sprintf(fnom,"A:%lx.OUT",time(NULL));
 #elif PC==2
  sprintf(fnom,"%lx OUT",time(NULL));
 #endif
 if((pf=fopen(fnom,"wb, recfm=v, lrecl=32756"))==NULL)
    {printf("NO PUC OBRIR EL FITXER DE SORTIDA"); exit(0);}
}
 
presentacio()
{char frase[60], data[30]; time_t lt, time(); short i;
 sprintf(frase,"MODEL D'APRENENTATGE HISTORIC-GEOGRAFIC");
 printf("%10s***  %s  ***\n"," ",frase);
 printf("Data:"); lt=time(NULL); sprintf(data,ctime(&lt)); printf(data);
 for(i=0;fnom[i];i++) fnom[i]=toupper(fnom[i]); imprimir(fnom);
 prch(' '); data[strlen(data)-1]=37;
 imprimir(data); imprimir(frase); prch(37);
}
 
imprimir(fila)
char *fila;
{short int L;
 for(L=0;L<strlen(fila);L++) prch(fila[L]);}
 
inicialitzar()
{short int N,U; char inici[20];
 MXz=(1<<mx)-1;
 printf("Factor d'adjustament=\n"); scanf("%d",&ka);
 sprintf(inici,"  ka=%d",ka); imprimir(inici);
 printf("Factor de resignacio=\n"); scanf("%d",&kr);
 sprintf(inici,"  kr=%d",kr); imprimir(inici);
 printf("Factor de naixement=\n"); scanf("%d",&kn);
 sprintf(inici,"  kn=%d",kn); imprimir(inici);
 printf("Factor d'impacte=\n"); scanf("%f",&ki);
 sprintf(inici,"  ki=%.2f",ki); imprimir(inici);
 prch(37);
 for(N=0;N<NP;N++)
           {m[N]=mx; Mz[N]=MXz; bz[N]=0; B[N]=0; neonat[N]=1;
            natal[N]=50+kn*(NP/2-abs(NP/2-abs(N-NP/4)));
            for(U=0;U<=MM;U++) P[U][N]=0;
            for(U=0;U<=Mz[N];U++) B[N]+=(F[U][N]=natal[N]*(U==U0));
           }
}
 
valors()
{short int U,j,Uj,N; float cost,forza;
 printf("\n");
 for(U=0;U<=MM;U++)
    {PI[U]=0; cost=1;
     for(j=0;j<mm;j++)
        {Uj=((U>>j)%2);
         PI[U]+=Uj;
         if(!Uj) {cost*=mm-1-j; cost/=mm-1;}
        }
     PI[U]/=mm;
     ferocitat[U]=1-cost; forza=(float)U/((1<<(mm-1))-1);
     st[U]=forza*ferocitat[U];
     em[U]=(float)U/MM;
    }
 /*
 printf("\nU \t PI[U] \t ferocitat[U] \t st[U] \t em[U]\n");
 for(U=0;U<=MM;U++)
    {printf("%d \t %.3f \t\t %.3f ",U,PI[U],ferocitat[U]);
     printf("\t %.3f \t %.3f\n",st[U],em[U]);}
 */
 for(N=0;N<NP;N++) repressio(N);
}
 
estadistica()
{int mj,Mj,U; double sqrt();
 for(mj=0;mj<=mm;mj++)
    {Mj=(1<<mj)-1; PR0[mj]=0; SR0[mj]=0;
     for(U=0;U<=Mj;U++) {PR0[mj]+=PI[U]; SR0[mj]+=PI[U]*PI[U];}
     PR0[mj]/=Mj+1; SR0[mj]/=Mj+1;
     SR0[mj]=sqrt((double)SR0[mj]-PR0[mj]*PR0[mj]);
    }
 /*
 printf("\nm \t PR0[m] \t SR0[m]\n");
 for(mj=0;mj<=mm;mj++)
  printf("%d \t %.3f \t %.3f\n",mj,PR0[mj],SR0[mj]);
 */
}
 
repressio(N)
short int N;
{short int U;
 for(U=0;U<=MM;U++) sts[U][N]=st[U];
}
 
pre_impr()
{/* inicialitzar impressora amb avanc automatic de linia, alta qualitat,
 paper continu de 12 polzades, interlinia 1/2, lletra estreta i subindex */
 prch(39); prch('@'); prch(39); prch('m'); prch(1); prch(39); prch('c');
 prch(39); prch('C'); prch(0); prch(12);
 prch(39); prch('A'); prch(46); prch(61);
 prch(39); prch('S'); prch(1); prch(39);
}
 
interrupcions()
{printf("\nInterrupcio en t=\n"); scanf("%ld",&Ti);}
 
pre_pant()
{short int i;
 printf("\nN. de Sistemes amb P(u/n) compres en els \n");
 printf("intervals:     0 ");
 for(i=0;i<mx;i++) printf("   a   %3d",100/(1<<(mx-i)));
 printf("   a   100");
}
 
iteracio()
{avaluacio();
 pantalla();
 impressora();
 if(T+100>Ti) sortida();
 T+=100;
 cientifi();
 aprenentatge();
 iteracio();
}
 
avaluacio()
{short int N,NN,U,d;
 for(N=0;N<NP;N++)
    if (B[N]!=0)
        {a[N]=.09; bz[N]=0; aT[N]=0; Tr[N]=0;
         for(U=0;U<=MXz;U++)
            {P[U][N]=(float)F[U][N]/B[N]; a[N]-=P[U][N]*U*.18/MM;
             aT[N]+=U*P[U][N]; Tr[N]+=(ferocitat[U]+1)*P[U][N];
             if(U>Mz[N]) bz[N]+=P[U][N];
            }
         for(NN=0;NN<NP;NN++)
            {d=NP/2-abs(NP/2-abs(N-NN)); IMPA[N][NN]=0;
             for(U=0;U<=MXz;U++) IMPA[N][NN]+=P[U][N]*imp[U][d]; }
         if(m[N]<mm) incremen(N);
         aT[N]/=ka*100*MM; Tr[N]*=kr*100;
        }
}
 
incremen(N)
short int N;
{if(rand()+32767*bz[N]+32767*B[N]/progres>=32767)
    {++m[N]; Mz[N]=(1<<m[N])-1; repressio(N);
     if(m[N]>mx) {mx=m[N]; MXz=Mz[N]; imprmx=1;} }
}
 
pantalla()
{short int RP[mm+1],j,U,N,R;
 if(imprmx==1) pre_pant();
 printf("\nt=%ld\n",T);
 for(U=0;U<=MXz;U++)
    {RP[0]=NP; for(R=mx;R>0;R--) RP[R]=0;
     for(N=0;N<NP;N++)
        {if(B[N]==0) --RP[0];
         else for(R=mx;R>0;R--)
                  if(P[U][N]*(1<<R)>=1) {++RP[mx+1-R]; --RP[mx-R];}
        }
     printf("estat U="); for(j=mx-1;j>=0;j--) printf("%d",(U>>j)%2);
     for(R=0;R<=mx;R++) printf("%10d",RP[R]); printf("\n");
    }
}
 
impressora()
{short int N,U; char buff[30],c[2];
 if(imprmx==1)
     {sprintf(buff,"dimensio maxima=%d",mx);
      imprimir(buff); prch(37); imprmx=0;}
 sprintf(buff,"t=%5ld    ",T); imprimir(buff);
 for(N=0;N<NP;N++)
    {if(B[N]!=0)
         {c[0]='-';
          for(U=0;U<=MXz;U++)
              if(P[U][N]>.5) {sprintf(c,"%X",U); break;}
         }
     else c[0]=' ';
     prch(c[0]);
    }
 if(T%1000==0)
    {prch(37); sprintf(buff,"   dim(N)= "); imprimir(buff);
     for(N=0;N<NP;N++)
         {if(B[N]>0) sprintf(c,"%X",m[N]); else c[0]=' '; prch(c[0]);}
    }
 /*  prch(39); prch('J'); prch(46); */ prch(37);
}
 
aprenentatge()
{short int U,V,N,NN,MXzz,S; double sqrt();
 float sigma,cosigma,PG[MM+1],PL[MM+1],PGM,SPG,L,sigm[MM+1],kk;
 K1=K2=0;
 for(N=0;N<NP;N++)
    {if(B[N]!=0)
     {L=0; PGM=0; SPG=0;
      for(U=0;U<=MXz;U++)
        {/* REPRESSIO */
         sigma=0; S=0;
         for(NN=0;NN<NP;NN++) if(B[NN]!=0)
            {++S;
             for(V=0;V<=MXz;V++) if(V!=U)
                {sigma+=P[V][NN]*P[V][NN]*sts[V][NN]*IMPA[NN][N];
                 if(sigma<0) {acuseta(sigma,NN,N); sortida();}
                }
            }
         if(S!=0) sigma/=S; else sortida();
/*       if(sigma>1) {acuseta(sigma,47,N); sortida();}
*/       cosigma=1-sigma; sigm[U]=sigma;
         sts[U][N]+=100*(sigma-sts[U][N])*aT[N];
         PG[U]=PI[U]*cosigma;
         /* COMUNICACIO CIENTIFICA */
         PL[U]=P[U][N];
         for(NN=0;NN<NP;NN++) if(NN!=N)
              if(B[NN]!=0) PL[U]+=EMz[N]*P[U][NN]*EMz[NN]*IMPA[NN][N];
         L+=PL[U]; PGM+=PG[U]*PL[U]; SPG+=PG[U]*PG[U]*PL[U];
         K1+=PI[U]*P[U][N]; K2+=sts[U][N]*P[U][N];
        }
      if(L<.999) {printf("\nATENCIO! L=%.3f<1\n",L); acuseta(sigma,N,NM);
                  sortida();}
      PGM/=L; SPG/=L; SPG-=PGM*PGM;
      if(SPG<0) {printf("\nATENCIO! SPG<0\n"); sortida();}
      SPG=sqrt((double)SPG);
      B[N]=0;
      /* RESIGNACIO */
      if(neonat[N])
         {PR[N]=PR0[m[N]]; SR[N]=SR0[m[N]];
          neonat[N]=0;}
      else
         {PR[N]+=100*(PGM-PR[N])/Tr[N];
          SR[N]+=100*(SPG-SR[N])/Tr[N];}
      /* APRENENTATGE INDIVIDUAL */
      for(U=0;U<=MXz;U++)
         {if(SR[N]>.0001)
          F[U][N]+=200*(PG[U]-PR[N])*PL[U]*SR0[m[N]]/SR[N];
          F[U][N]*=(F[U][N]>0);
          B[N]+=F[U][N];
         }
      /* REVERSIBILITAT DE LA DESTRUCCIO ECOLOGICA */
      if(B[N]==0)
        {if(S==1) {printf("\nHOLOCAUSTO!\n"); sortida();}
         kk=(float)S/--S;
         if(kk*k<1) {for(U=0;U<=MXz;U++) PI[U]*=kk; k*=kk;}
         else {for(U=0;U<=MXz;U++) PI[U]/=k; k=1;}
        }
      NM=N;
     }
     /* REPRODUCCIO */
     if(rand()+32767*a[N]+32767*(1-a[N])*B[N]/tanatos>=32767)
        {if(B[N]==0)
           {short NN=N; while(B[NN]==0) NN=(NN+1)%NP;
            m[N]=m[NN]; Mz[N]=(1<<m[N])-1;}
         for(U=0;U<=MXz;U++) F[U][N]=natal[N]*(U<=Mz[N]);
         B[N]=natal[N]*(Mz[N]+1);
         repressio(N);
         neonat[N]=1;
        }
    }
    /* acuseta(sigma,NM,50); */
    MXzz=(MXz<MM/2)?MXz:MM/2;
    printf("\nPG  =");
    for(U=MXz/2;U<=MXz;U++) printf("%.3f ",PG[U]);
    printf("\nsts =");
    for(U=MXz/2;U<=MXz;U++) printf("%.3f ",sts[U][NM]);
    printf("\nN=%d m=%d U=%d a %d  ",NM,m[NM],MXz/2,MXz);
    printf("PGM=%.3f SPG=%.3f  PR=%.3f SR=%.3f\n",PGM,SPG,PR[NM],SR[NM]);
    /* DESTRUCCIO ECOLOGICA */
    if(K1+K2>NP)
       {kk=K1/(NP-K2); for (U=1;U<=MXz;U++) PI[U]/=kk; k/=kk;
        printf("          K=%.0f+%.0f=%.0f",K1,K2,K1+K2);
       }
    if(k<1) printf("          ecologia=%.2f",k);
}
 
cientifi()
{short N,U;
 for(N=0;N<NP;N++) if(B[N]!=0)
    {EMz[N]=0;
     for(U=0;U<=MXz;U++) EMz[N]+=P[U][N]*em[U];}
}
 
sortida()
{void exit(); char control=0; long i;
 printf("\nPer terminar, pulsa la tecla 't' i RETURN\n");
 printf("Per continuar, pulsa la tecla 'c' i RETURN\n");
 while(control!='t' && control!='c') control=getchar();
 if(control=='t')
    {acuseta(0.,NM,NM);
     printf("\n***  %s  ***\n\n",fnom);
     /* prch(39); prch('@'); prch(39); prch('J'); prch('l'); */ exit(0);}
 interrupcions();
 for(i=0;i<100000;i++);
}
 
prch(c)
char c;
{fputc(c,pf);}
 
impacte()
{int U,d,dm; float forza,mimp,Mimp;
 for(d=0;d<=NP/2;d++)
 {for(U=1;U<=MM;U++)
    {forza=(float)2*U/(MM-1);
     mimp=1-ki+ki*(forza>1)*(forza-1)*(MM-1)/(MM+1);
     Mimp=2-mimp;
     dm=NP/2; if(forza<1) dm*=forza;
     imp[U][d]=mimp+(dm>d)*(Mimp-mimp)*(dm-d)/dm;
    }
  imp[0][d]=2*(!d);
 }
}
 
acuseta(sigma,N,NN)
float sigma; short N,NN;
{short mj=m[N],U;
 printf("\n  sigma=%f k=%.2f",sigma,k);
 printf("\n  IMPA[%d][%d]=%.3f",N,NN,IMPA[N][NN]);
 printf("\nm=%d U=",mj); for(U=mj-1;U<=MXz;U++) printf("%4d ",U);
 printf("\n   PI="); for(U=mj-1;U<=MXz;U++) printf("%.2f ",PI[U]);
 printf("\n    P="); for(U=mj-1;U<=MXz;U++) printf("%.2f ",P[U][N]);
 printf("\n  sts="); for(U=mj-1;U<=MXz;U++) printf("%.2f ",sts[U][N]);
 printf("\n");
}
 
mostrar()
{int d,U,i;
 printf("\n\n\n");  for(i=0;i<10000;i++);
 for(d=0;d<=NP/2;d++) {printf("\n%2d: ",d);
     for(U=0;U<8;U++) printf("%.3f ",imp[U][d]);}
 printf("\n");
 for(d=0;d<=NP/2;d++) {printf("\n%2d: ",d);
     for(U=8;U<16;U++) printf("%.3f ",imp[U][d]);}
}