/* $Header: /disks/hermosa/public/NEURON/CA1-99/bp/lib/TEST/newshiftsyn.c 2000/02/25 poirazi */
/* Adding a temporal offset for stimulation of trains */
/* $Header: /disks/newport/freeway/brannon/rs/p-and-s/SourceCode/C/newshiftsyn.c,v 1.7 1999/03/24 18:19:21 brannon Exp $ */
/* $Log: newshiftsyn.c,v $
* Revision 1.7 1999/03/24 18:19:21 brannon
* Ok, now NEURON should be happy.
*
* Revision 1.6 1999/03/24 18:09:25 brannon
* The end of newshiftsyn files is not compatible with NEURON's vector
* scanf() function. Working on this.
*
* Revision 1.5 1998/10/22 10:05:03 brannon
* The compiler croaked on nested comments so I fixed my attempt to
* simply comment out and copy the old code by removing the old code and
* associated comments.
*
* Revision 1.4 1998/10/22 10:02:58 brannon
* I changed the sprintf which created filenames to write the p_seed in
* integer as opposed to float from.
*
* Revision 1.3 1998/05/07 19:33:44 brannon
* rm shiftsyn.log -> rm -f shiftsyn.log
* My runs were held up because it was waiting for a reply to removal of
* a previous shiftsyn.log
*
* Revision 1.2 1998/04/08 18:33:29 brannon
* Write a \n after each dt of synaptic input.
*
* Revision 1.1 1998/04/08 18:32:48 brannon
* Initial revision
*
* Revision 1.12 1997/04/15 19:35:59 niebur
* *** empty log message ***
*
* Revision 1.11 1995/02/01 01:55:28 ernst
* Eliminated finite loop condition: if s=1, the restart in the main loop after
* enforce_refrac failed (cond=1) was ineffective, because make_train_offset
* and make-spike trains each time gave identical results. This happened because
* the sigmas in them contained terms 1-s and thus became zero.
*
* This has been corrected by restarting the spiketrain now COMPLETELY
* (ie, all synapses) when that occurs.
*
* PS this condition occurred for " shiftsyn testfile 100 500 0.1 100 1 0.2 6074"
*
* Revision 1.10 1995/01/28 22:00:46 ernst
* Writing now a log file at the beginning of the run which is
* deleted at its successful completion.
*
* Revision 1.9 1995/01/28 08:04:27 ernst
* Now enforcing refractory period (in new routine enforce_refrac).
* Also found out that the Numerical recipes generator has a problem
* if its seed is too large (> 54000, apparently). The reason is
* unknown, but I do not allow such large seeds anymore.
*
* Revision 1.8 1995/01/20 22:13:26 ernst
* Added more debugging message in context with the consistency check in write_to_disk.
*
* Revision 1.7 1995/01/18 02:41:29 ernst
* Added a global variable ("errfile") which is identical to outfilename and
* which is printed out each time there is an error message. This way, it will
* be possible to determine which run makes a problem in a series of many runs
* in which the shell might garble the order of messages.
*
* Revision 1.6 1995/01/18 02:01:03 ernst
* removed debugging messages
*
* Revision 1.5 1994/10/29 02:25:07 ernst
* Fixed a problem in make_spike_trains: the sigma of the Gaussian there
* was dependend on the total time. There is no reason this should be
* the case.
*
* Revision 1.4 1993/10/04 08:11:56 ernst
* Introduced computation and writing to a file of
* the total activity at a given time (variable totact).
*
* Revision 1.3 1993/09/23 01:17:08 ernst
* Introduced larger scatter for shift between spike trains.
* The scatter is now in units of n_tot, NOT in units of the
* average spiking frequency.
* */
/* $Id: newshiftsyn.c,v 1.7 1999/03/24 18:19:21 brannon Exp $ */
/***************************************************************************/
/* */
/* Purpose of this program: make time-dependent input for use with NEURON */
/* */
/* Input: explained before 'main' */
/* */
/* Output: file: */
/* */
/* syndata<parameters> contains spike times for the given number */
/* of synapses (read by NEURON), */
/* */
/* Note: the random number generator is defined in /cit/ernst/num-rec.h */
/* */
/****************************************************************************/
#include "/usr/include/stdlib.h"
#include "/usr/include/string.h"
//#include "/usr/include/math.h"
#include "ken.h"
#include "num-rec.h"
#define RAND_MAX 2147483647
#define myrand ( (double) random() / RAND_MAX)
char errfile[1000]; /* global var for debugging purposes*/
/********************************************************************/
/* */
/* Allocate array space for times (doubles) */
/* */
/********************************************************************/
double *alloctimes(int nelements)
{
int i;
double *p;
p = calloc(nelements,sizeof(double));
if(!p)
{
printf("allocation failure in allocspikes, working on %s\n; aborting",
errfile);
exit(1);
}
for(i=0;i<nelements;i++)
p[i]=0;
return(p);
}
/********************************************************************/
/* */
/* Allocate array space for spikes (ints) */
/* */
/********************************************************************/
int *allocspikes(int nelements)
{
int i;
int *p;
p = calloc(nelements,sizeof(int));
if(!p)
{
printf("allocation failure in allocspikes, working on %s\n; aborting",
errfile);
exit(1);
}
for(i=0;i<nelements;i++)
p[i]=0;
return(p);
}
/********************************************************************/
/* */
/* Read input from console */
/* */
/********************************************************************/
void read_input(int argc,char *argv[],int *p_n_syn,double *p_t_tot,double *p_dt,double *p_av_rate,double *p_p,double *p_s,int *p_seed, double *p_offset, char *outfilename,char *totactname)
{
float f_t_tot,f_dt,f_av_rate,f_p,f_s,f_offset;
if(argc != 10)
{
printf
("Usage: shiftsyn outfile n_syn t_tot dt av_rate s p seed offset\n (Error file is %s)\n",errfile);
exit(1);
}
else
{
sscanf(argv[2],"%i",p_n_syn);
sscanf(argv[3],"%g",&f_t_tot);
sscanf(argv[4],"%g",&f_dt);
sscanf(argv[5],"%g",&f_av_rate);
sscanf(argv[6],"%g",&f_s);
sscanf(argv[7],"%g",&f_p);
sscanf(argv[8],"%i",p_seed);
*p_offset = atof(argv[9]);
}
/* (translating the floats into doubles)*/
if(f_p<0 || f_p>1 || f_p<0 || f_s>1 )
{
printf("p and s must be from [0..1];, working on %s EXIT!\n\n",errfile);
exit(1);
}
*p_t_tot = f_t_tot;
*p_dt = f_dt ;
*p_p=f_p;
*p_s=f_s;
/*
Brannon: I dont have the faintest idea why p_seed would be cast as float here.
But I do have a non-faint idea as to why I want to make it integer.
*/
/* making the complete outfile name (NB: dt is in ms):*/
sprintf(outfilename,"%s-%d-%.2f-%.2f-%.2f-%.2f-%.2f-%d-%.2f",\
argv[1],*p_n_syn,f_t_tot,f_dt,f_av_rate,f_s,f_p,*p_seed,*p_offset);
/* also write the complete outfile name to global variable which is
used in error messages as a signature:*/
sprintf(errfile,"%s-%d-%.2f-%.2f-%.2f-%.2f-%.2f-%d-%.2f",\
argv[1],*p_n_syn,f_t_tot,f_dt,f_av_rate,f_s,f_p,*p_seed,*p_offset);
/* making the name for the total activity file:*/
sprintf(totactname,"%s-%d-%.2f-%.2f-%.2f-%.2f-%.2f-%d-%.2f",\
"totact",*p_n_syn,f_t_tot,f_dt,f_av_rate,f_s,f_p,*p_seed,*p_offset);
/* t_tot and dt are in ms, therefore also the rate has to be in spikes */
/* per ms and not spikes per second: */
*p_av_rate = f_av_rate/1000.;
}
/********************************************************************/
/* */
/* make_train_offset: */
/* Compute the offset for the whole spiketrain. */
/* A Gaussian random number with sigma s/(1-s) */
/* */
/* s=0 <==> all spiketrains are completely periodic, */
/* s=1 <==> all spiketrains are completely aperiodic, */
/* s=0.5 <==> spiketrains are distributed with sigma= */
/* */
/********************************************************************/
double make_train_offset(double p, double s, double av_rate, \
double t_tot,int *p_seed)
{
double offset,delta_t;
double sigma_p;
delta_t = 1./av_rate;
if(s > APRECISION && s<=1.) /* finite value for s */
{
sigma_p = p*(1-s)/s * delta_t; /* determine sigma of gaussian */
offset= sigma_p*gasdev(p_seed);
}
else if(fabs(s) < APRECISION) /* s is zero ==> "infinite sigma" */
{
offset= p* t_tot * ran1(p_seed); /* pick random offset */
}
else /* error! */
nrerror("error in make_train_offset: s must be 0 <= s <= 1");
return(offset);
}
/********************************************************************/
/* */
/* dsort: */
/* sorting routine, from numrec. */
/* Identical to routine sort there, exc. for one difference (shown) */
/* */
/********************************************************************/
void dsort(n,ra)
int n;
double ra[];
/* declaration of ra and rra is the only difference to numrec's sort */
{
int l,j,ir,i;
double rra;
l=(n >> 1)+1;
ir=n;
for (;;) {
if (l > 1)
rra=ra[--l];
else {
rra=ra[ir];
ra[ir]=ra[1];
if (--ir == 1) {
ra[1]=rra;
return;
}
}
i=l;
j=l << 1;
while (j <= ir) {
if (j < ir && ra[j] < ra[j+1]) ++j;
if (rra < ra[j]) {
ra[i]=ra[j];
j += (i=j);
}
else j=ir+1;
}
ra[i]=rra;
}
}
/********************************************************************/
/* */
/* sort_train: */
/* bring spikes in range [0..t_tot] and */
/* sort spiketrain in ascending order. */
/* */
/********************************************************************/
sort_train(double *spiketrain,int length_train,double t_tot)
{
int i;
for (i=0;i<length_train; i++)
{
/*bring in range [0 .. t_tot]:*/
spiketrain[i]=fmod(spiketrain[i],t_tot);
while(spiketrain[i]<0)
spiketrain[i] += t_tot;
if(spiketrain[i]>t_tot || spiketrain[i] <0)
{
printf("sort_train: This cannot happen\n, working on %s",errfile);
exit(1);
}
}
/* DEBUG:*/
/*
printf("\nNext train, in range 0-%g, unordered (working on %s):\n",
t_tot,errfile);
for (i=0;i<length_train; i++)
{
printf("%i %g\n",i,(float)spiketrain[i]);
}
*/
/* use numerical recipes heapsort routine. ATTENTION: subtract 1 from
the pointer to the array spiketrain because the routine expects array
indices 1 to n, not 0 to n-1. See NRC p. 15 (unit offsets vs. zero
offsets) */
dsort(length_train,spiketrain-1);
}
/********************************************************************
* *
* enforce_refrac: *
* enforce refractory period. *
* *
* Method: If a spike is closer than t_ref to the next one, *
* the next one is pushed by t_refrac. This MAY lead the last *
* spike(s) to be pushed beyond t_tot. If that happens, we push *
* back the last spikes by t_ref/2 (not t_ref; this guarantees *
* that the order remains conserved). *
* *
* returns 1 if successful, -1 if not *
* *
********************************************************************/
int enforce_refrac(double *spiketrain, int syn, int length_train,
double t_ref, double dt,double t_tot)
{
int i,j,done;
if(t_ref <= dt) {
printf("ERROR1 in enforce_refrac, working on %s\n; aborting\n",
errfile);
exit(1);
}
done=0;
while( !done) { /* loop goes to l_t-2; see below for very last interval*/
for (i=0;i<length_train-2; i++)
if(spiketrain[i+1]-spiketrain[i] <= t_ref) { /* if too close */
spiketrain[i+1] += t_ref; /* push next spike*/
j=i+1; /*if spike order violated, keep pushing till ok:*/
if(j+1 >= length_train)
return(-1); /* failure; make new spike train*/
while((spiketrain[j]) > spiketrain[j+1]) {
spiketrain[j+1] += t_ref;
j++;
if(j+1 >= length_train) /*make sure j-loop stays w/in array bounds*/
return(-1); /*else: failure; make new spike train*/
}
}
/*Now spikes are separated by at least t_ref. Check if sth. fell off end:*/
done=1;
for (i=1;i<length_train; i++)
if(spiketrain[i] > t_tot) { /*if yes, put the last one back a bit:*/
spiketrain[i-1] -= t_ref/2.; /*subtract only t_r/2 to keep ordered*/
if(spiketrain[i-1]<0) /*failure; make new spike train*/
return(-1);
done=0; /*have to do the while loop again (possibly more than once)*/
}
}
/* The only interval that has not been checked is the very last one.
In the unlikely case this is a problem, we would have to do a
major re-ordering (spikes are already ordered!).
Therefore, if this happens, declare failure, return with (-1) and make
an entirely new spike train:*/
if(spiketrain[length_train-1]-spiketrain[length_train-2] <= t_ref)
return(-1);
else {
/*Debug:*/
/*
printf("\nOrdered train %i, in range 0-%g (working on %s):\n",
syn,t_tot,errfile);
for (i=0;i<length_train; i++)
{
printf("%i %g\n",i,(float)spiketrain[i]);
}
*/
return(1);
}
}
/********************************************************************/
/* */
/* append_train: */
/* append spiketrain to alltrains. */
/* In this array, spiketrains are arranged one after the other. */
/* ie., first the complete train for synapse 1, then for syn. 2 etc */
/* */
/********************************************************************/
append_train(double *sortedtrain,double *alltrains, int syn,\
int length_train)
{
int i;
for (i=0;i<length_train; i++)
{
/*store all synapses together in array:*/
alltrains[syn*length_train+i]=sortedtrain[i];
}
}
/********************************************************************/
/* */
/* write_to_disk: */
/* write spike trains to disk and */
/* write the total activity to a separate file. */
/* also do consistency checking. */
/* */
/********************************************************************/
void write_to_disk(double *alltrains,FILE *outfile,FILE *f_totact,\
int length_train,double t_tot,int n_syn,double dt,double offset)
{
double t;
int nspike,syn,outsyn;
int *indnexttime;
double *nexttime;
int totact;
FILE *problem_file;
char consist_filename[1000];
/*To avoid that we have to read the array for each synapse
completely at each time step, we store the next time at which a
spike occurs in nexttime and its index in the array in
indnexttime:*/
indnexttime = allocspikes(n_syn);
nexttime = alloctimes(n_syn);
for(syn=0;syn<n_syn;syn++) /* find first spike for each synapse */
{
indnexttime[syn]=0;
nexttime[syn]=alltrains[syn*length_train]+offset; /* time of first spike */
}
for(t=0;t<=t_tot; t+=dt)
{
totact=0; /*total activity at this time*/
for (syn=0;syn<n_syn;syn++)
{
if(fabs(t-nexttime[syn])<dt) /*this synapse has a spike at time t */
{
fprintf(outfile,"1"); /* write spike to disk */
indnexttime[syn] ++; /* next spike of this synapse: */
nexttime[syn] = offset + alltrains[syn*length_train+indnexttime[syn]];
totact ++; /*one spike more in total activity*/
if(nexttime[syn] > t_tot)
{
nexttime[syn] = t_tot+dt;
}
}
else /* synapse has no spike at time t */
{
fprintf(outfile,"0");
}
if (syn<(n_syn-1)) {
fprintf(outfile," ");
}
}
if (t < (t_tot-dt)) {
fprintf(outfile,"\n");
}
/*write total activity for this time:*/
// fprintf(f_totact,"%g %d\n",t,totact);
}
/*make sure all spikes have been found (consistency check):*/
// for(syn=0;syn<n_syn;syn++)
// {
// if(indnexttime[syn] != length_train)
// {
// printf("Working on %s: Someth. wrong w. synapse %d:\n",errfile,syn);
// printf("Last spike: %i != length of spiketrain: %i\n",
// indnexttime[syn],length_train);
// sprintf(consist_filename,"ERROR_%s",errfile);
// problem_file = fopen(consist_filename,"a");
//
// fprintf(problem_file,"\n\n*******\n\nn_syn: %i\n",n_syn);
// fprintf(problem_file,"length_train: %i\n",length_train);
// fprintf(problem_file,"t_tot: %i\n",t_tot);
// fprintf(problem_file,"syn: %i\n\n",syn);
// for (outsyn=0;outsyn<n_syn;outsyn++)
// for(nspike=0;nspike<length_train;nspike++)
// fprintf(problem_file,"%i %i: %g\n",
// outsyn,nspike,alltrains[outsyn*length_train+nspike]);
// fclose (problem_file);
// }
// }
}
/********************************************************************/
/* */
/* make_spike_train: */
/* add shift to individual spikes */
/* given by gaussian with sigma_s */
/* */
/********************************************************************/
void make_spike_trains(double *spiketrain,double *gentrain,double p,\
double s,double av_rate,double t_tot,double toff,\
int *p_seed)
{
double t,delta_t;
double sigma_s;
int i;
double alpha;
alpha=2.;
delta_t = 1./av_rate;
i=0;
/*sigma_s = t_tot; this is wrong, no reason for t_tot dependence!*/
sigma_s = (1-p)*(1-s)*alpha*delta_t; /* determine sigma of gaussian */
for(t=0.;t<t_tot;t+=delta_t) /* make spike train with gaussian */
{
spiketrain[i] = gentrain[i] + toff + sigma_s*gasdev(p_seed);
i++;
}
}
/********************************************************************/
/* */
/* make_generator_train: */
/* make master spike train */
/* */
/********************************************************************/
void make_generator_train(double *gentrain,double t_tot,\
double av_rate,double p,int *p_seed)
{
double t,delta_t;
double sigma_g;
int i;
delta_t = 1./av_rate;
i=0;
if(p >APRECISION && p<=1.) /* finite value for p */
{
sigma_g = (1-p)/p * delta_t; /* determine sigma of gaussian */
for(t=0.;t<t_tot;t+=delta_t) /* make spike train with gaussian */
{
gentrain[i] = t+sigma_g*gasdev(p_seed);
i++;
}
}
else if(fabs(p) < APRECISION) /* p is zero ==> "infinite sigma" */
{
for(t=0.;t<t_tot;t+=delta_t)
{
gentrain[i] = t_tot * ran1(p_seed); /* pick random spike times */
i++;
}
}
else /* error! */
nrerror("error in make_generator_train: p must be 0 <= p <= 1");
}
/********************************************************************/
/* */
/* Main. */
/* Parameters: */
/* n_syn - Number of synapses */
/* t_tot - total simulation time (milliseconds) */
/* dt - simulation time step (milliseconds) */
/* av_rate - average spike rate (spikes per second) */
/* p - periodicity parameter, 0<p<1 */
/* s - synchronicity parameter, 0<s<1 */
/* offset - spike temporal offset (milliseconds) */
/* */
/********************************************************************/
main(int argc, char *argv[])
{
int n_syn,syn;
int length_train;
int start_offset;
double t_tot,dt,av_rate, offset;
FILE *outfile,*f_totact,*logfile;
char outfilename[1000],totactname[1000];
double *spiketrain, *gentrain, *alltrains;
double p,s;
double toff;
int seed,cond;
int i;
double t_ref=1.; /*refractory period, in ms*/
strcpy(errfile,"shiftsyn.err");
read_input(argc, argv,&n_syn,&t_tot,&dt,&av_rate,&p,&s,&seed,&offset,outfilename,totactname);
// printf( "%g\n", offset);
OPENFILE(logfile,"shiftsyn.log","w",exit(1));
fprintf(logfile, "Right now, shiftsyn is working on %s\n",errfile);
fclose(logfile);
if(seed>50000) {
printf("seed too large; problem for Num. Recipes random generator\n");
exit(1);
}
if(t_tot*av_rate<=1)
nrerror("run too short (less than one spike on average), EXIT");
else
length_train = floor(t_tot*av_rate);
spiketrain = alloctimes(length_train);
gentrain = alloctimes(length_train);
alltrains = alloctimes(n_syn*length_train);
outfile = fopen(outfilename, "w");
// f_totact = fopen(totactname, "w");
make_generator_train(gentrain,t_tot,av_rate,p,&seed);
for(syn=0;syn<n_syn;syn++) /* loop over all synapses */
{
/* compute global offset for whole spiketrain of this synapse: */
toff = make_train_offset(p,s,av_rate,t_tot,&seed);
/* add shifts to individual spike trains */
make_spike_trains(spiketrain,gentrain,p,s,av_rate,t_tot,toff,&seed);
/* sort spiketrain */
sort_train(spiketrain,length_train,t_tot);
/* enforce refractory period */
cond=enforce_refrac(spiketrain,syn,length_train,t_ref,dt,t_tot);
if (cond == 1) { /*normal case*/
/* append to alltrains:*/
append_train(spiketrain,alltrains,syn,length_train);
}
else if (cond == -1) { /*make a new spiketrain for same syn*/
if(s<0.9) {
syn --; /*if s!=1 it is enough to restart this synapse */
}
else { /*if s==1 need also new gen. train and all new synapses */
make_generator_train(gentrain,t_tot,av_rate,p,&seed);
syn=-1;
}
}
else {
printf("ERROR in main: this cannot happen\n");
exit(1);
}
}
/* All spiketimes are now computed. Write them to disk:*/
write_to_disk(alltrains,outfile,f_totact,length_train,t_tot,n_syn,dt,offset);
/* Debugging output: */
/*
for(i=0;i<n_syn*length_train; i++)
{
printf("Working on %s: %d %f\n",errfile,i,alltrains[i]);
printf("%f \n",alltrains[i]);
}
*/
/* successfully completed, rm logfile: */
system("rm -f shiftsyn.log");
}