#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "MT.h"
#include <unistd.h>
#define getrandom(max1) ((genrand_int32()%(max1)))
#define N 1000
#define M 100
#define D 20
int Ne = 800;
int Ni;
float sm = 10.0;
int post[N][M];
float s[N][M], sd[N][M];
short delays_length[N][D];
short delays[N][D][M];
int N_pre[N], I_pre[N][3*M], D_pre[N][3*M];
float *s_pre[N][3*M], *sd_pre[N][3*M];
float LTP[N][1001+D], LTD[N];
float a[N], c[N], d[N];
float v[N], u[N];
int N_firings;
const int N_firings_max = 100*N;
int firings[N_firings_max][2];
double Ap = 0.1;
double Am = 0.12;
double tau = 20.0;
double window;
void format(double *a1, double *a2, double *c1, double *c2, double *d1, double *d2, int *time, int *Ne, int *seed, double *Ap, double *Am, double *tau, int argc, char **argv)
{
int i;
for(i=0; i<argc/2; i++){
if(argv[2*i+1] == NULL) break;
else if (argv[2*(i+1)] == NULL){
fprintf(stderr, "param format [-a1,-a2,-c1,-c2,-d1,-d2,-time,-Ne,-seed,-Ap,-Am,-tau] [num]\n");
exit(1);
}
else if(strcmp(argv[2*i+1], "-a1") == 0){
*a1 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-a2") == 0){
*a2 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-c1") == 0){
*c1 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-c2") == 0){
*c2 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-d1") == 0){
*d1 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-d2") == 0){
*d2 = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-time") == 0){
*time = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-Ne") == 0){
*Ne = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-seed") == 0){
*seed = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-Ap") == 0){
*Ap = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-Am") == 0){
*Am = atof(argv[2*(i+1)]);
}
else if(strcmp(argv[2*i+1], "-tau") == 0){
*tau = atof(argv[2*(i+1)]);
}
else {
fprintf(stderr, "param format [-a1,-a2,-c1,-c2,-d1,-d2,-time,-Ne,-seed,-Ap,-Am,-tau] [num]\n");
exit(1);
}
}
}
void initialize(double a1, double a2, double c1, double c2, double d1, double d2){
int i,j,k,jj,dd, exists, r;
Ni = N - Ne;
window = exp(-1/tau);
for (i=0;i<Ne;i++){
a[i] = a1;
c[i] = c1;
d[i] = d1;
}
for (i=Ne;i<N;i++){
a[i] = a2;
c[i] = c2;
d[i] = d2;
}
for (i=0;i<N;i++){
for (j=0;j<M;j++){
do{
exists = 0;
if (i<Ne){
r = getrandom(N);
}
else{
r = getrandom(Ne);
}
if (r==i){
exists=1;
}
for (k=0;k<j;k++){
if (post[i][k]==r){
exists = 1;
}
}
}
while (exists == 1);
post[i][j]=r;
}
}
for (i=0;i<Ne;i++){
for (j=0;j<M;j++){
s[i][j] = 6.0;
}
}
for (i=Ne;i<N;i++){
for (j=0;j<M;j++){
s[i][j]=-5.0;
}
}
for (i=0;i<N;i++){
for (j=0;j<M;j++){
sd[i][j]=0.0;
}
}
for (i=0;i<N;i++){
short ind = 0;
if (i<Ne){
for (j=0;j<D;j++){
delays_length[i][j] = M/D;
for (k=0; k<delays_length[i][j]; k++)
delays[i][j][k] = ind++;
}
}
else{
for (j=0;j<D;j++){
delays_length[i][j] = 0;
}
delays_length[i][0] = M;
for (k=0; k<delays_length[i][0]; k++){
delays[i][0][k] = ind++;
}
}
}
for (i=0;i<N;i++){
N_pre[i] = 0;
for (j=0; j<Ne; j++){
for (k=0; k<M; k++){
if (post[j][k] == i){
I_pre[i][N_pre[i]]=j;
for (dd=0;dd<D;dd++){
for (jj=0;jj<delays_length[j][dd];jj++){
if (post[j][delays[j][dd][jj]]==i){
D_pre[i][N_pre[i]]=dd;
}
}
}
s_pre[i][N_pre[i]]=&s[j][k];
sd_pre[i][N_pre[i]++]=&sd[j][k];
}
}
}
}
for (i=0;i<N;i++){
for (j=0;j<1+D;j++){
LTP[i][j]=0.0;
}
}
for (i=0;i<N;i++){
LTD[i]=0.0;
}
for (i=0;i<N;i++){
v[i] = -65;
}
for (i=0;i<N;i++){
u[i] = 0.2*v[i];
}
N_firings = 1;
firings[0][0] = -D;
firings[0][1] = 0;
}
int main(int argc, char **argv)
{
int i, j, k, sec, t;
int seed = getpid();
int time = 3600;
float I[N];
FILE *fs;
double a1 = 0.02;
double a2 = 0.1;
double c1 = -65.0;
double c2 = -65.0;
double d1 = 8.0;
double d2 = 2.0;
format(&a1, &a2, &c1, &c2, &d1, &d2, &time, &Ne, &seed, &Ap, &Am, &tau, argc, argv);
init_genrand(seed);
initialize(a1, a2, c1, c2, d1, d2);
for (sec=0; sec<time; sec++){
for (t=0; t<1000; t++){
for (i=0; i<N; i++){
I[i] = 0.0;
}
for (k=0; k<N/1000; k++){
I[getrandom(N)]=20.0;
}
for (i=0; i<N; i++){
if (v[i]>=30){
v[i] = c[i];
u[i] += d[i];
LTP[i][t+D] = Ap;
LTD[i] = Am;
for (j=0; j<N_pre[i]; j++){
*sd_pre[i][j] += LTP[I_pre[i][j]][t+D-D_pre[i][j]-1];
}
firings[N_firings ][0] = t;
firings[N_firings++][1] = i;
if (N_firings == N_firings_max) {
N_firings=1;
}
}
}
k = N_firings;
while (t-firings[--k][0] < D){
for (j = 0; j < delays_length[firings[k][1]][t-firings[k][0]]; j++){
i = post[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]];
I[i] += s[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]];
if (firings[k][1] <Ne){
sd[firings[k][1]][delays[firings[k][1]][t-firings[k][0]][j]]-=LTD[i];
}
}
}
for (i=0;i<N;i++){
v[i]+=0.5*((0.04*v[i]+5)*v[i]+140-u[i]+I[i]);
v[i]+=0.5*((0.04*v[i]+5)*v[i]+140-u[i]+I[i]);
u[i]+=a[i]*(0.2*v[i]-u[i]);
LTP[i][t+D+1] = window*LTP[i][t+D];
LTD[i] *= window;
}
}
for (i=1;i<N_firings;i++){
if (firings[i][0] >= 0){
printf("%d %d\n", firings[i][0] + 1000*sec, firings[i][1]);
}
}
for (i=0;i<N;i++){
for (j=0;j<D+1;j++){
LTP[i][j]=LTP[i][1000+j];
}
}
k=N_firings-1;
while (1000-firings[k][0]<D){
k--;
}
for (i=1; i<N_firings-k; i++){
firings[i][0] = firings[k+i][0]-1000;
firings[i][1] = firings[k+i][1];
}
N_firings = N_firings-k;
for (i=0;i<Ne;i++){
for (j=0;j<M;j++){
s[i][j]+=0.01+sd[i][j];
sd[i][j]*=0.9;
if (s[i][j]>sm){
s[i][j]=sm;
}
if (s[i][j]<0){
s[i][j]=0.0;
}
}
}
}
return 0;
}