/* Program do generacji polaczen synaptycznych miedzy subnets
Na wyjsciu ma produkowac plik z polaczeniami z zapisana para neuronow waga i
delay. Nastepny program (lub druga czesc tego samego) moze wykorzystac to
do generacji struktury neuronow i synaps we wlasciwym formacie.
01/29/03 P. Kudela, losuj2 prawdopodobienstwo polaczenia zanika exp.
a w zapisz2 delay wzrasta liniowo
05/29/03 P. Kudela, now multi-class supported
*/
#define RISC
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>
#include <math.h>
#include "lnet.h"
#define randomm(num) (int)((double)num*(double)random()/((double)(RAND_MAX+1.0)))
#define randomize() srandom((unsigned)time(NULL))
int lrodz; /* number of neuron classes */
int *NN; /* NN[lrodz] total number of neurons within each class */
int *NOF; /* NOF[lrodz] table of index offsets */
int **N; /* N[lrodz][3] size of subnetwork for each class */
int **k; /* k[lrodz][3] scaling factor to align class networks */
int ***o; /* o[lrodz][lrodz][3] neighborhood array */
int **s; /* s[lrodz][lrodz] synapsy na we/wy z otoczenie */
int **w,**sdw; /* w[lrodz][lrodz] sdw[][] synaptic weights and dispersions */
int **d,**sdd; /* d[lrodz][lrodz] sdd[][] synaptic delays and dispersions */
double **Ap,**lambda;/* Ap[lrodz][lrodz] lambda[][] amplitude and spatial decay rate of connection probability */
int *tab; /* temporary neighborhood table tab[maxot]*/
double *dist, *costheta, *sintheta; /* table of distances between candidate neurons and costheta, sintheta, updated 08/23/05 */
int *postx, *posty; /* Array of postsynaptic neuron coordinates for accepted target */
double *zclass; /* Array of class z-values within cortical model */
int *horaxonflagpos,*veraxonflagpos, *axonflag; /* Array of flags for stimulated axon portions */
int *horaxonflagneg,*veraxonflagneg;
/* Neuron Class Type Firing Pattern
1 Layer II/III Pyramid RS
2 Layer IV Spiny RS
3 Layer V Pyramid IB
4 Layer VI Pyramid IB
5 Basket Layers II-VI FS
6 Double Bouquet II-VI FS
7 Chandelier FS */
#define P23 1
#define P4 2
#define P5 3
#define P6 4
#define BASKET 5
#define BOUQUET 6
#define CHANDELIER 7
#define BASKETANG 0.0 //1.5708 /* for offset costheta distr for basket cell axons */
#define PIIIANG 0.0 /* for offset costheta distr for PIII cell axons */
#define PVANG 0.0 /*for offset costheta distr for PV cell axons */
#define PVIANG 0.0 /* for offset costheta distr for PVI cell axons */
#define SPATIALSCALE 25 /* center to center column spacing in microns */
#define BASKETLAMBDA 50 /* 1/e decay for basket cell axon connections */
#define BASKETNEAR 50 /* distance in microns within which basket conn'ions are isotropic */
#define PIIINEARRAD 300 /* Near connection radius for PIII neurons */
#define STRIPWIDTH 500 /* 500 micron strip width for domain simulations */
#define PIIIFARAVE 2000 /* Ave dist to long range connection spot */
#define PIIISPOTRAD 200 /* Long range spot radius */
#define BOUNDARYLAYER 300 /* Boundary layer distance, to reduce connections */
#define ISLENGTH 30 /* One half initial segment length (microns) */
/* Definitions for E-field effects */
#define SEG 500 /* spacing btwn nodes of Ranvier for 5.7 micron myelinated axon */
//#define VOLT .2 /* stimulation voltage applied to metal electrode in voltage controlled mode */
#define I0 .010 /* Stimulation current in current controlled mode in Amps */
#define SIGMA 0.3 /* Tissue conductivity in A/Vm */
#define THRESHD2V 2 /* Threshold for activation fct (mV), second difference of voltage along fiber element */
#define A 1000 /* Electrode disk radius (microns) */
#define STIMPOS 2340 /* Electrode Z-value (microns), taken with z=0 at the top of the white matter */
#define DISKX 32 /* Center of electrode disk x-value (lattice constants) */
#define DISKY 32 /* Center of electrode disk y-value (lattice constants) */
#define PI 3.14159
/* Cortical layer z-values */
#define ZVIB 162 /* VIbeta */
#define ZVIA 488 /* VIalpha */
#define ZV 1130 /* V */
#define ZIVCBETA 1230 /* IVCbeta */
#define ZIVCALPHA 1230 /* IVCalpha */
#define ZIVB 1230 /* IVB */
#define ZIVA 1230 /* IVA */
#define ZIII 1540 /* III */
#define ZII 1940 /* II */
#define ZI 2185 /* I */
#define Gi 7.6e-07 /* Intracompartment conductivity in inverse Ohms */
struct SYN_TMP
{
NEURON_INDEX neuron_pre;
NEURON_INDEX neuron_post;
NEURON_INDEX preclass_num;
WEIGHT weight;
DELAY delay;
VAR axdelay;
};
int sort_delay(const void *, const void *);
int sort_neu_num(const void *, const void *);
void set_neu(struct NEURON *);
unsigned seed_save();
int main(int argc, char *argv[])
{
int i,j,l,jp;
int nr,ir,ix,iy,iz; /* index for neurons*/
int maxot, tmpot;
int xyz[3]={0,0,0};
int no,nl;
ulong ns,is;
int *neu_pre;
int Nneu;
int torus_flag;
int syn_n;
char label[]="xyz",name[10];
FILE *inputfile,*outputfile,*graphfile,*axflagfile;
struct SYN_TMP *syn_tmp;
struct SYNAPSE syn;
struct NEURON *neu;
int otoczenie(int,int,int*,int);
int otoczenie2(int,int,int*,int);
int losuj(int*,int);
int losuj2(int*,double*,double*,double*,int,double,double,int,int,int*);
extern int zapisz(int,int,int,int,FILE *);
extern int zapisz2(int,int,int,int,double,FILE *,FILE *,FILE *,int);
extern int sort_neu_num(const void *, const void *);
int nsub;
int pre[3];
double horzmid,horxmid,horymid,horrhomid,horVmid,horrhopre,horVpre,horrhopost,horVpost;
double horseglength,hordelsquaredV;
double verzmid,verxmid,verymid,verrhomid,verVmid,verrhopre,verVpre,verrhopost,verVpost;
double verseglength,verdelsquaredV,aradius,coeff;
int synapsesum;
double currentinjsum, avecurrentinj;
int nps,s_j,s_i,critdex;
/*
1.Reading and checking input file information
*/
/* printf("Made it to file checking\n"); */
if(argc<3)
{
fprintf(stderr,"USAGE:\n\tmk_link <file_name> <sqrt(no of subnets in all)> [<seed>]\n");
return 1;
}
if(argc>3)
srandom((unsigned)atoi(argv[3]));
else
printf("seed %u\n",seed_save()); /* seed for srand */
srand48(127L);
//inputfile = fopen(strcpy(name,argv[1]),"r");
inputfile = fopen(strcat(strcpy(name,argv[1]),".net"),"r");
if(inputfile==NULL){fprintf(stderr, "Network definition file not present %s\n",name);return -1;}
fscanf(inputfile,"%d",&torus_flag);
fscanf(inputfile,"%d",&lrodz);
nsub=atoi(argv[2]);
/* Memory allocations */
N=(int**)malloc(lrodz*sizeof(int*)); /* table allocation for size of networks for each class */
for(i=0;i<lrodz;i++)
N[i]=(int*)malloc(3*sizeof(int));
k=(int**)malloc(lrodz*3*sizeof(int)); /* table allocation for network scaling factor for each class */
for(i=0;i<lrodz;i++)
k[i]=(int*)malloc(3*sizeof(int));
o=(int ***)malloc(lrodz*sizeof(int **)); /* allocating neighborhood tables */
for(i=0;i<lrodz;i++)
o[i]=(int **)malloc(lrodz*sizeof(int *));
for(i=0;i<lrodz;i++)
for(j=0;j<lrodz;j++)
o[i][j]=(int *)malloc(3*sizeof(int));
s=(int **)malloc(lrodz*sizeof(int*)); /* synaptic table allocation */
for(i=0;i<lrodz;i++)
s[i]=(int*)malloc(lrodz*sizeof(int));
w=(int **)malloc(lrodz*sizeof(int*)); /* allocation for synaptic weights */
for(i=0;i<lrodz;i++)
w[i]=(int*)malloc(lrodz*sizeof(int));
sdw=(int **)malloc(lrodz*sizeof(int*)); /* allocation for synaptic weight dispersion */
for(i=0;i<lrodz;i++)
sdw[i]=(int*)malloc(lrodz*sizeof(int));
d=(int **)malloc(lrodz*sizeof(int*)); /* allocation for synaptic delay */
for(i=0;i<lrodz;i++)
d[i]=(int*)malloc(lrodz*sizeof(int));
sdd=(int **)malloc(lrodz*sizeof(int*)); /* allocation for synaptic delay dispersion */
for(i=0;i<lrodz;i++)
sdd[i]=(int*)malloc(lrodz*sizeof(int));
Ap=(double **)malloc(lrodz*sizeof(double)); /* used to compute probability of connection */
for(i=0;i<lrodz;i++)
Ap[i]=(double*)malloc(lrodz*sizeof(double));
lambda=(double **)malloc(lrodz*sizeof(double)); /* used to commpute probability of connection */
for(i=0;i<lrodz;i++)
lambda[i]=(double*)malloc(lrodz*sizeof(double));
zclass=(double*)malloc(lrodz*sizeof(double));
/* Load cortical layer z-information for E-field calc, array index by class number */
zclass[0]= ZII;
zclass[1]= ZIVB;
zclass[2]= ZV;
zclass[3]= ZVIA;
zclass[4]= ZII;
zclass[5]= ZII;
zclass[6]= ZIVB;
/* initialize mapping and scaling arrays */
for(i=0;i<lrodz;i++)
{
N[i][0]=N[i][1]=nsub; //square mapping
N[i][2]=1;
k[i][0]=k[i][1]=k[i][2]=1;
}
/* input parameters from file to fill arrays */
for(i=0;i<lrodz;i++)
for(j=0;j<lrodz;j++)
{
fscanf(inputfile,"%d,%d,%d",&o[j][i][0],&o[j][i][1],&o[j][i][2]);
fscanf(inputfile,"%d",&s[j][i]);
fscanf(inputfile,"%d",&w[j][i]);
fscanf(inputfile,"%d",&sdw[j][i]);
fscanf(inputfile,"%d",&d[j][i]);
fscanf(inputfile,"%d",&sdd[j][i]);
fscanf(inputfile,"%lf",&Ap[j][i]);
fscanf(inputfile,"%lf",&lambda[j][i]);
}
fclose(inputfile);
/*
2. Offset calculations and table filling
*/
/* printf("Made it to offset calculations\n"); */
NN=(int*)malloc(lrodz*sizeof(long int)); /* table with number of subnetworks containing each class */
for(i=0;i<lrodz;i++)
NN[i]=(N[i][0]?N[i][0]:1)*(N[i][1]?N[i][1]:1)*(N[i][2]?N[i][2]:1);
NOF=(int*)malloc(lrodz*sizeof(int)); /* table with offsets */
for(i=1,NOF[0]=0;i<lrodz;i++)
NOF[i]=NOF[i-1]+NN[i-1];
/* maximal neigborhood calculation */
maxot=0;
if(torus_flag!=2) /* i.e. not choosing from whole array */
for(i=0;i<lrodz;i++)
for(j=0;j<lrodz;j++)
{
tmpot=1;
for(l=0;l<3;l++)
tmpot*=2*o[i][j][l]+1;
if(tmpot > maxot)
maxot=tmpot;
}
else
for(i=0;i<lrodz;i++)
{
if(maxot<NN[i])
maxot=NN[i];
}
tab=(int*)malloc(maxot*sizeof(int));
dist=(double*)malloc(maxot*sizeof(double));
costheta=(double*)malloc(maxot*sizeof(double));
sintheta=(double*)malloc(maxot*sizeof(double));
postx=(int*)malloc(maxot*sizeof(int));
posty=(int*)malloc(maxot*sizeof(int));
horaxonflagpos=(int*)malloc(maxot*sizeof(int));
veraxonflagpos=(int*)malloc(maxot*sizeof(int));
horaxonflagneg=(int*)malloc(maxot*sizeof(int));
veraxonflagneg=(int*)malloc(maxot*sizeof(int));
axonflag=(int*)malloc(maxot*sizeof(int));
Nneu=0;
for(i=0;i<lrodz;i++)
Nneu+=NN[i]; /* calculate subnetworks X # of classes */
for(i=0;i<lrodz;i++)
printf("\nNN[%d]=%d, NOF[%d]=%d",i,NN[i],i,NOF[i]); /* co by sprawdzic */
printf("\nmaxot=%d,Nneu=%d\n",maxot,Nneu);
/*
3.Setting up for output of links
*/
outputfile=fopen("syn.tmp","wb");
graphfile=fopen("syn.graph","wb");
axflagfile=fopen("syn.axflag","wb");
if(outputfile==NULL){fprintf(stderr,"\nError in opening synapse output file");exit(1);}
if(graphfile==NULL){fprintf(stderr,"\nError in opening synapse output file");exit(1);}
if(axflagfile==NULL){fprintf(stderr,"\nError in opening axflagfile file");exit(1);}
/*randomize();*/
nr=0; /* subnetwork "neuron" number */
ns=0; /* synapse index */
currentinjsum=0; /* in microamps */
synapsesum=0;
for(ir=0;ir<lrodz;ir++) /* loop on classes of pre-neurons */
{
//printf("Working on Preclass %d\n",ir+1);
for (ix=0,xyz[0]=0;ix<N[ir][0];ix++,xyz[0]+=k[ir][0]) /* x-loop */
{
//printf("ix= %d\n",ix+1);
for(iy=0,xyz[1]=0;iy<N[ir][1];iy++,xyz[1]+=k[ir][1]) /* y-loop */
for(iz=0,xyz[2]=0;iz<N[ir][2];iz++,xyz[2]+=k[ir][2],nr++) /* z-loop */
for(jp=0;jp<lrodz;jp++) /* loop on classes of post-neurons */
{
/* generate tables with numbers of neurons of each class jp,in the neighborhood of neuron nr (neuron from class ir
connects to jp (used to be the other way around) */
switch(torus_flag)
{
case 0: /* no torus boundary condition */
no=otoczenie(nr,ir,xyz,jp);
break;
case 1: /* torus boundary */
no=otoczenie2(nr,ir,xyz,jp);
break;
case 2: /* case 2 choosing from whole array */
for(i=NOF[jp],no=0;i<NN[jp]+NOF[jp];i++)
if(i!=nr)
tab[no++]=i;
break;
default:
fprintf(stderr,"wrong option for torus flag %d\n",torus_flag);
exit(0);
}
if(!no)continue;
l=(no<s[jp][ir])?no:s[jp][ir];
/* if no is greater than the number of allowed synapses, then set l=# of allowed synapses */
//printf("Allowed/required conn comp no=%d,preclass=%d,postclass=%d,s[jp][ir]=%d\n",no,ir+1,jp+1,s[jp][ir]);/* network linkage checks */
/* Boundary effect correction, reducing allowed connections for neurons in specific column number ranges, and row number ranges */
nps=nsub; //sqrt(subnets/node)*sqrt(nodes)
s_j=(nr)%nps; /* column index of presynaptic subnetwork */
s_i=(nr)/nps-ir*nps; /* row index of presynaptic subnetwork */
critdex=BOUNDARYLAYER/SPATIALSCALE;
//printf("Subnet indices: column=%d row=%d nps=%d critdex=%d nr=%d\n",s_j,s_i,nps,critdex,nr);
//printf("Pre-l: l=%d\n",l);
if ( (s_j<critdex) && (s_i>=critdex) && (s_i<=((nps-1)-critdex)) ){
l=(l*5)/8;
}
if ( (s_j>((nps-1)-critdex)) && (s_i>=critdex) && (s_i<=((nps-1)-critdex)) ){
l=(l*5)/8;
}
if ( (s_i<critdex) && (s_j>=critdex) && (s_j<=((nps-1)-critdex)) ){
l=(l*5)/8;
}
if ( (s_i>((nps-1)-critdex)) && (s_j>=critdex) && (s_j<=((nps-1)-critdex)) ){
l=(l*5)/8;
}
if ( (s_j<critdex) && (s_i<critdex) ){
l=(l*3)/8;
}
if ( (s_j<critdex) && (s_i>((nps-1)-critdex)) ){
l=(l*3)/8;
}
if ( (s_j>((nps-1)-critdex)) && (s_i<critdex) ){
l=(l*3)/8;
}
if ( (s_j>((nps-1)-critdex)) && (s_i>((nps-1)-critdex)) ){
l=(l*3)/8;
}
//printf("Post-l: l=%d\n",l);
/* Form connections and determine stimulation locations */
for(i=0;i<l;i++) /* losuje z calego otoczenia - nowa wersja */
{ /* z losuj2 01/29/03 */
for(j=0;j<3;j++)
pre[j]=xyz[j]/k[jp][j];
// printf("no=%d,preclass=%d,postclass=%d,l=%d,i=%d\n",no,ir+1,jp+1,l,i); /* network linkage checks */
nl=-1;
while (nl==-1){
nl=losuj2(tab,dist,costheta,sintheta,no,Ap[jp][ir],lambda[jp][ir],jp,ir,pre);
}
//printf("i= %d,l= %d,nl= %d\n",i,l,nl);
if(nl!=-1)
{
/* Set synapse flag if gradient requirement met */
/* Check E field gradient requirement */
/* Charged balance, bipolar pulse */
/* Horizontal axon segment */
aradius=A;
aradius=aradius/(1000*1000);
coeff=I0/(4*PI*SIGMA*aradius);
//printf("\n aradius=%f, coeff=%f",aradius,coeff);
axonflag[nl]=0;
horaxonflagpos[nl]=0;
horaxonflagneg[nl]=0;
horzmid=zclass[jp];
horxmid=(pre[0]+postx[nl])/2;
horymid=(pre[1]+posty[nl])/2;
horrhomid=sqrt(pow((DISKX-horxmid),2)+pow((DISKY-horymid),2))*25;
horVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-horzmid),2)+pow((horrhomid-A),2))+sqrt(pow((STIMPOS-horzmid),2)+pow((horrhomid+A),2))));
horrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
horVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopre-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopre+A),2))));
horrhopost=sqrt(pow((DISKX-postx[nl]),2)+pow((DISKY-posty[nl]),2))*25;
horVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopost-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopost+A),2))));
horseglength=sqrt(pow((postx[nl]-pre[0]),2)+pow((posty[nl]-pre[1]),2))*25;
/* hordelsquaredV in mV */
hordelsquaredV = (horVpre+horVpost-2*horVmid)*1000;
if (hordelsquaredV > THRESHD2V)
horaxonflagpos[nl]=1;
if (-hordelsquaredV > THRESHD2V)
horaxonflagneg[nl]=1;
// Vertical axon segment (Old Version)
//veraxonflagpos[nl]=0;
//veraxonflagneg[nl]=0;
//verzmid=(zclass[ir]+zclass[jp])/2;
//verxmid=pre[0];
//verymid=pre[1];
//verrhomid=sqrt(pow((DISKX-verxmid),2)+pow((DISKY-verymid),2))*25;
//verVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid-A),2))+sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid+A),2))));
//verrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
//verVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre-A),2))+sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre+A),2))));
//verrhopost=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
//verVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((verrhopost-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((verrhopost+A),2))));
//verseglength=sqrt(pow((zclass[jp]-zclass[ir]),2));
///* verdelsquaredV in mV */
//verdelsquaredV = (verVpre+verVpost-2*verVmid)*1000;
//if (verdelsquaredV > THRESHD2V)
// veraxonflagpos[nl]=1;
//if (-verdelsquaredV > THRESHD2V)
// veraxonflagneg[nl]=1;
// Vertical axon segment (IS Version)
veraxonflagpos[nl]=0;
veraxonflagneg[nl]=0;
verzmid=(2*zclass[ir]-2*ISLENGTH)/2;
verxmid=pre[0];
verymid=pre[1];
verrhomid=sqrt(pow((DISKX-verxmid),2)+pow((DISKY-verymid),2))*25;
verVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid-A),2))+sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid+A),2))));
verrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
verVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre-A),2))+sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre+A),2))));
verrhopost=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
verVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]+2*ISLENGTH),2)+pow((verrhopost-A),2))+sqrt(pow((STIMPOS-zclass[ir]+2*ISLENGTH),2)+pow((verrhopost+A),2))));
verseglength=2*ISLENGTH; //sqrt(pow((zclass[jp]-zclass[ir]),2));
/* verdelsquaredV in mV */
verdelsquaredV = (verVpre+verVpost-2*verVmid)*1000;
if (verdelsquaredV > THRESHD2V)
veraxonflagpos[nl]=1;
if (-verdelsquaredV > THRESHD2V)
veraxonflagneg[nl]=1;
/* if (-verdelsquaredV > THRESHD2V){
currentinjsum+=(-verdelsquaredV)*Gi;
synapsesum=synapsesum+1;
} */
//if ( (horaxonflagpos[nl]==1) || (veraxonflagpos[nl]==1) )
if ( veraxonflagpos[nl]==1 )
axonflag[nl]=1;
//if ( (horaxonflagneg[nl]==1) || (veraxonflagneg[nl]==1) )
if ( veraxonflagneg[nl]==1 )
axonflag[nl]=2;
ns+=zapisz2(nr,tab[nl],jp,ir,dist[nl],outputfile,graphfile,axflagfile,axonflag[nl]);
/* printf("\npostindex=%d, postclass=%d, preindex=%d, preclass=%d",nr,ir+1,tab[nl],jp+1);*/ /* network linkage checks */
// tab[nl]=-1;
}
} /* koniec petli losowania */
} /* jp */
}
}
fclose(outputfile);
fclose(graphfile);
fclose(axflagfile);
// avecurrentinj=currentinjsum/synapsesum;
//printf("currentinjsum=%f,avecurrentinj=%f,synapsesum=%i\n",currentinjsum,avecurrentinj,synapsesum); /* network linkage checks */
printf("\n");
return 0;
} // koniec main()
//sort_delay - uzywana w qsort() przy porzadkowaniu delay'ow
int sort_delay(const void * a, const void * b )
{
char delay1[1],delay2[1];
struct SYN_TMP *syn1=NULL, *syn2=NULL;
syn1=(struct SYN_TMP *)a;
syn2=(struct SYN_TMP *)b;
*delay1=syn1->delay;
*delay2=syn2->delay;
return( strcmp(delay1,delay2) );
}
/*
sort_neu_num - uzywana w qsort()
*/
int sort_neu_num(const void * a, const void * b)
{
NEURON_INDEX neu1,neu2;
struct SYN_TMP *syn_tmp1=NULL,*syn_tmp2=NULL;
syn_tmp1=(struct SYN_TMP *)a;
syn_tmp2=(struct SYN_TMP *)b;
neu1=syn_tmp1->neuron_pre;
neu2=syn_tmp2->neuron_pre;
return(neu1-neu2);
}
/*
4. Szczegoly generacji tablicy otoczenia:
*/
int rozrzut(int sd)
{
int los;
los=randomm(100);
if(los < 50 )
return( -1*randomm(sd) );
else
return( randomm(sd) );
}
int losuj(int *tablica,int index)
{
int i,j;
do
{
i=randomm(index);
}while(tablica[i]<0);
j=tablica[i];
tablica[i]=-1;
return j;
}
int losuj2(int *tablica, double *distance, double *costht,double *sintht,int index, double Ap, double lambda, int classpost, int classpre,int *precoord)
/* zwraca index w tab, a nie numer neuronu jak losuj() */
{
int i,j,flagspatial;
//double Ap=1.0;
//double lambda=6;
double los,losang,los1,los2,cutoff;
double p3nearrad,p3farave,p3spotrad;
double stripsize,preyy,postycoord,prefloor,postfloor;
double precheck,postcheck;
// i=randomm(index);
do
{
i=randomm(index);
// printf("i=%d,tablica[i]=%d,index=%d\n",i,tablica[i],index); /* network linkage checks */
}while(tablica[i]<0);
switch(classpre)
{
case 0: /* Pre - Layer II/III Pyramid */
/* Patch Handling */
p3nearrad=PIIINEARRAD/SPATIALSCALE;
stripsize=STRIPWIDTH/SPATIALSCALE;
//p3farave=PIIIFARAVE/SPATIALSCALE;
//p3spotrad=PIIISPOTRAD/SPATIALSCALE;
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
//losang=drand48();
// los2=1;
flagspatial=0;
if (distance[i]<=p3nearrad)
flagspatial=1;
else {
/* Stripflag A=even, B=odd */
preyy=precoord[1];
postycoord=posty[i];
prefloor=floor(preyy/stripsize);
postfloor=floor(postycoord/stripsize);
precheck=prefloor-2*floor(prefloor/2.0);
postcheck=postfloor-2*floor(postfloor/2.0);
//printf("precheck=%f,postcheck=%f,preyy=%f\n",precheck,postcheck,preyy);
if ( precheck == postcheck)
flagspatial=0;
else
flagspatial=0;
//printf("precheck=%f,postcheck=%f,flagspatial=%d\n",precheck,postcheck,flagspatial);
}
/* Angular isotropy routine */
//los1=(costht[i]*cos(PIIIANG)+sintht[i]*sin(PIIIANG));
//los1*=los1;
//los1*=los1;
//los1*=los1;
/*printf("p3nearrad=%lf,p3farave=%lf,p3spotrad=%lf,flagspatial=%i,distance=%lf\n",p3nearrad,p3farave,p3spotrad,flagspatial,distance[i]);*/
if( ( flagspatial ) /*&& (losang < los1)*/ )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
//los=drand48();
//// los2= Ap*exp(-1.0*distance[i]/lambda);
//los2=1;
////printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
//if( los < los2 )
//{
// return i;
//}
//else
//{
// // tablica[i]=-1;
// return -1;
//}
//break;
case 1: /* Pre - Layer IV Spiny */
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
los2=1;
//printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
if( los < los2 )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
case 2: /* Pre - Layer V Pyramid */
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
//losang=drand48();
los2=1;
//los1=(costht[i]*cos(PVANG)+sintht[i]*sin(PVANG));
//los1*=los1;
//los1*=los1;
//los1*=los1;
//printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
if( ( los < los2 ) /* && (losang < los1) */ )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
case 3: /* Pre - Layer VI Pyramid */
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
//losang=drand48();
los2=1;
//los1=(costht[i]*cos(PVIANG)+sintht[i]*sin(PVIANG));
//los1*=los1;
//los1*=los1;
//los1*=los1;
//printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
if( ( los < los2 ) /* && (losang < los1) */ )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
case 4: /* Pre - Basket Cell */
cutoff = BASKETNEAR/SPATIALSCALE;
if ( distance[i] <= cutoff ){ /* Isotropic connections to 50 microns */
return i;
}
else {
los=drand48();
//losang=drand48();
lambda = BASKETLAMBDA/SPATIALSCALE;
los2= Ap*exp(-1.0*distance[i]/lambda);
//los1=(costht[i]*cos(BASKETANG)+sintht[i]*sin(BASKETANG));
//los1*=los1;
//los1*=los1;
//los1*=los1;
// los2=1;
/* printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2); */
if( (los < los2) /* && (losang < los1) */ )
{
/* printf("index=%i, neuron=%i, costht=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,costht[i],losang,los1);*/
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
}
break;
case 5: /* Pre - Double Bouquet */
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
los2=1;
//printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
if( los < los2 )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
case 6: /* Pre - Chandelier Cell */
los=drand48();
// los2= Ap*exp(-1.0*distance[i]/lambda);
los2=1;
//printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);
if( los < los2 )
{
return i;
}
else
{
// tablica[i]=-1;
return -1;
}
break;
default:
fprintf(stderr,"Problem with neuron class in losuj2 %d\n",classpre);
exit(0);
}
}
int zapisz(int inr, int inl, int i, int j, FILE * outputfile)
{
struct SYN_TMP syn_tmp;
extern int rozrzut(int);
syn_tmp.neuron_post=inl;
syn_tmp.neuron_pre=inr;
syn_tmp.weight=w[i][j]+rozrzut(sdw[i][j]);
syn_tmp.delay =d[i][j]+rozrzut(sdd[i][j]);
//printf("%i %i %i %i\n",syn_tmp.neuron_pre, syn_tmp.neuron_post , syn_tmp.weight,syn_tmp.delay );
return(fwrite(&syn_tmp,sizeof(struct SYN_TMP),1,outputfile));
}
int zapisz2(int inr, int inl, int i, int j, double distance, FILE * outputfile, FILE * graphfile, FILE * axflagfile, int axflag)
{
struct SYN_TMP syn_tmp;
extern int rozrzut(int);
int axonal_del=5; /* was 5 */
syn_tmp.neuron_post=inl;
syn_tmp.neuron_pre=inr;
syn_tmp.preclass_num=j;
syn_tmp.weight=w[i][j]+rozrzut(sdw[i][j]);
syn_tmp.delay =d[i][j]+rozrzut(sdd[i][j]);
//printf("distance=%3.2lf, axonal del=%3.2lf, delay=%i, delay + axonal del=%i\n", distance, distance*axonal_del, syn_tmp.delay, syn_tmp.delay+(int)rint(distance*axonal_del));
syn_tmp.axdelay=(int)rint(distance*axonal_del)+3*rozrzut(sdd[i][j]);
//printf("%i %i %i %i %i\n",syn_tmp.neuron_pre, syn_tmp.neuron_post , syn_tmp.weight, syn_tmp.delay , syn_tmp.axdelay);
//printf("syn_tmp.axdelay=%i\n", syn_tmp.axdelay);
if ((j)==4){
fwrite(&inr,sizeof(ushort),1,graphfile);
fwrite(&j,sizeof(ushort),1,graphfile);
fwrite(&inl,sizeof(ushort),1,graphfile);
fwrite(&i,sizeof(ushort),1,graphfile);
}
/* printf("axflag= %i\n",axflag); */
fwrite(&axflag,sizeof(ushort),1,axflagfile);
return(fwrite(&syn_tmp,sizeof(struct SYN_TMP),1,outputfile));
}
int otoczenie(int nr ,int ir,int xyz[3],int jp)
{
int nx,ny,no,i[3],ix,iy,iz,mini[3],maxi[3], *ot;
int j;
/* tab must be declared earlier with sufficient size */
ot = o[jp][ir];
/* obliczamy najpierw wspolrzedne i[3] neuronu na siatce jp, ktory jest najblizej x,y,z */ /*???*/
no=0;
for(j=0;j<3;j++)
i[j]=xyz[j]/k[jp][j];
/* teraz otoczenie
wersja bez zawijania (brzegowe maja mniejsze otoczenie) */ /*???*/
for(j=0;j<3;j++)
{
mini[j]=i[j]-ot[j];
maxi[j]=i[j]+ot[j]+1;
if(mini[j]<0)mini[j]=0;
if(maxi[j]>N[jp][j])maxi[j]=N[jp][j];
}
for(ix=mini[0];ix<maxi[0];ix++)
{ nx=ix*N[jp][1]*N[jp][2]+NOF[jp];
for(iy=mini[1];iy<maxi[1];iy++)
{ ny=nx+iy*N[jp][2];
for(iz=mini[2];iz<maxi[2];iz++)
{
if((j=ny+iz)!=nr){
dist[no]=sqrt((ix-i[0])*(ix-i[0])+(iy-i[1])*(iy-i[1])+(iz-i[2])*(iz-i[2]));
costheta[no]=(ix-i[0])/dist[no];
sintheta[no]=(iy-i[1])/dist[no];
/* printf("ix=%i, iy=%i, i[0]=%i,i[1]=%i,costheta=%lf,sintheta=%lf,dist=%lf\n",ix,iy,i[0],i[1],costheta[no],sintheta[no],dist[no]);*/
// printf("%6.1lf,",dist[no]); /* tablica z odleglosciami - nowa wersja 01/29/03 */
postx[no]=ix;
posty[no]=iy;
tab[no++]=j;}
}
}
}
//printf("\n\n");
return no;
}
int otoczenie2(int nr, int ir,int xyz[3],int jp)
{
int nx,ny,no,i[3],ix,iy,iz,mini[3],maxi[3], *ot;
int j,ii,NX,NY,NZ;
/* tab must be declared earlier with sufficient size */
ot = o[jp][ir];
/* obliczamy najpierw wspolrzedne i[3] neuronu na siatce jp, ktory jest najblizej x,y,z */ /* ??? */
no=0;
for(j=0;j<3;j++)
i[j]=xyz[j]/k[jp][j];
/* Torus case */
for(j=0;j<3;j++)
{
mini[j]=i[j]-ot[j];
maxi[j]=i[j]+ot[j]+1;
}
NX=N[jp][0];
NY=N[jp][1];
NZ=N[jp][2];
for(ix=mini[0];ix<maxi[0];ix++)
{ ii=(ix+NX)%NX;
nx=ii*NZ*NY+NOF[jp];
for(iy=mini[1];iy<maxi[1];iy++)
{ ii=(iy+NY)%NY;
ny=nx+ii*NZ;
for(iz=mini[2];iz<maxi[2];iz++)
{
if((j=ny+(iz+NZ)%NZ)!=nr)
tab[no++]=j;
}
}
}
return no;
}
void set_neu(struct NEURON * neu)
{
int j,kk;
switch(neu->param)
{
default:
neu->flags|=HISTO;
neu->polarity = 0x0;
for (kk=0;kk<NEXCCLASS;kk++){
neu->Ve_o[kk]=0;
neu->Ve_d[kk]=0;
}
for (kk=0;kk<NINHCLASS;kk++){
neu->Vi_o[kk]=0;
neu->Vi_d[kk]=0;
}
neu->interval=0;
neu->sum_interval=0;
for(j=0;j<SPIKE_LIST_LENGTH;j++)
neu->spikes[j]=0;
neu->first=EMPTY_LIST;
neu->last=0x1;
neu->V=V_0;
neu->W=W_0;
#ifdef CALCIUM
for(j=1;j<7;j++)
{
neu->C[j]=0.;
neu->U[j]=0.;
}
// neu->C=C_0;
neu->X=X_0;
neu->B=B_0;
#endif
}
}
unsigned seed_save()
{
unsigned seed;
seed = (unsigned)time(NULL);
srandom(seed);
return seed;
}