#include <stdlib.h>
#include <string.h>
/*
n is the number of values to return
ncon is the number of connections each guy gets
index[n][ncon] is the list of indices to which it connects
weight[n][ncon] is the list of weights
root is the address of lowest entry in the variable array
To wit:
for sparse type
name(i) produces values[i]
value[i] = sum(k=0..ncon-1)of(w[i][k]*root[index[i][k]])
ODE FILE CALL:
special f=sparse( n, ncon, w,index,rootname)
special f=conv(type,n,ncon,w,rootname)
index,w are tables that are loaded by the "tabular"
command.
rootname here is the name of the root quantity.
type is either ep0
for conv.
conv0 conve convp type
value[i]= sum(j=-ncon;j<=ncon;i++){
k=i+j;
0 type if(k>=0 && k<n)
value[i]+=wgt[j+ncon]*rootname[k]
e type k=abs(i+j); if(k<2n)
if(k>=n)k=n-k;
...
p type
k=mod(2*n+i+j,n)
for example discretized diffusion
tabular dd % 3 -1 1 3*abs(t)-2
special diff=conv(even, 51, 2, dd,v0)
v[0..50]'=f(v[j],w[j])+d*diff([j])
another example
nnetwork
tabular wgt % 51 -25 25 .032*cos(.5*pi*t/25)
special stot=conv(0, 51, 25, wgt,v0)
v[0..50]'=f(v[j],w[j])-gsyn*stot([j])*(v([j])-vsyn)
last example -- random sparse network 51 cells 5 connections each
tabular w % 251 0 250 .2*rand(1)
tabular con % 251 0 250 flr(rand(1)*5)
special stot=sparse(51,5,w,con,v0)
v[0..50]'=f(v[j],w[j])-gsyn*stot([j])*(v([j])-vsyn)
more stuff:
special f=mmult(n,m,w,root) -- note that here m is the number of return
values and n is the size of input root
f[j] = sum(i=0,n-1)of(w(i+n*j)*root(i)) j=0,..,m-1
special f=fmmult(n,m,w,root1,root2,fname)
f[j] = sum(i=0,n-1)of(w(i+n*j)*fname(root1[i],root2[j])
special f=fsparse( n, ncon,w,index,root1,root2,fname)
special f=fconv(type,n,ncon,w,root1,root2,fname)
sum(j=-ncon,ncon)w(j)*fname(root1[i+j],root2[i])
similarly for fsparse
special k=fftcon(type,n,wgt,v0)
uses the fft to convolve v0 with the wgt which must be of
length n if type=periodic
2n if type=0
special f=interp(type,n,root)
this produces an interpolated value of x for
x \in [0,n)
type =0 for linear (all that works now)
root0,....rootn-1 are values at the integers
Like a table, but the integer values are variables
*/
#include <math.h>
#include <stdio.h>
extern int NODE;
#define MAX_TAB 50
#define IC 2
extern int fftn (int /* ndim */,
const int /* dims */[],
double /* Re */[],
double /* Im */[],
int /* isign */,
double /* scaling */);
/* simple network stuff */
typedef struct {
double xlo,xhi,dx;
double *y,*x;
int n,flag,interp,autoeval;
int xyvals;
/* flag=0 if virgin array, flag=1 if already allocated; flag=2 for function
interp=0 for normal interpolation, interp=1 for 'step'
table and finally, xyvals=1 if both x and y vals are needed (xyvals=0
is faster lookup )*/
char filename[128],name[12];
}TABULAR;
extern TABULAR my_table[MAX_TAB];
typedef struct {
int type,ncon,n;
char name[20];
int root,root2;
int f[20];
int iwgt;
int *gcom; /* for group commands */
double *values,*weight,*index;
double *fftr,*ffti,*dr,*di;
} NETWORK;
#define CONVE 0
#define CONV0 1
#define CONVP 2
#define SPARSE 3
#define FCONVE 10
#define FCONV0 11
#define FCONVP 12
#define FSPARSE 13
#define FFTCON0 4
#define FFTCONP 5
#define MMULT 6
#define FMMULT 7
#define GILLTYPE 25
#define INTERP 30
extern double variables[];
#define MAXNET 50
char *get_first(/* char *string,char *src */);
char *get_next(/* char *src */);
double evaluate();
NETWORK my_net[MAXNET];
int n_network=0;
double net_interp(double x, int i)
{
int jlo=(int)x;
double *y;
int n=my_net[i].n;
double dx=x-(double)jlo;
y=&variables[my_net[i].root];
if(jlo<0 || jlo>(n-1))return 0.0; /* out of range */
return (1-dx)*y[jlo]+dx*y[jlo+1];
}
double network_value(x, i)
double x;
int i;
{
int j=(int)x;
if(my_net[i].type==INTERP){
return net_interp(x,i);
}
if(j>=0&&j<my_net[i].n)
return my_net[i].values[j];
return 0.0;
}
init_net(double *v,int n)
{
int i;
for(i=0;i<n;i++)
v[i]=0.0;
}
add_spec_fun(name,rhs)
char *name;
char *rhs;
{
int i,ind,elen;
int type;
int iwgt,iind,ivar,ivar2;
int ntype,ntot,ncon,ntab;
char *str,cc;
char junk[256];
char rootname[20],wgtname[20],indname[20];
char root2name[20],fname[20];
type=is_network(rhs);
if(type==0)return 0;
printf("type=%d \n",type);
for(i=0;i<n_network;i++)
if(strcmp(name,my_net[i].name)==0)break;
ind=i;
if(ind>=n_network){
printf(" No such name %s ?? \n",name);
return 0;
}
switch(type){
case 1: /* convolution */
get_first(rhs,"(");
str=get_next(",");
ntype=-1;
if(str[0]=='E')ntype=CONVE;
if(str[0]=='0'||str[0]=='Z')ntype=CONV0;
if(str[0]=='P')ntype=CONVP;
if(ntype==-1){
printf(" No such convolution type %s \n",str);
return 0;
}
str=get_next(",");
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(")");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
my_net[ind].values=(double *)malloc((ntot+1)*sizeof(double));
init_net(my_net[ind].values,ntot);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].n=ntot;
my_net[ind].ncon=ncon;
printf(" Added net %s type %d len=%d x %d using %s var[%d] \n",
name,ntype,ntot,ncon,wgtname,ivar);
return 1;
break;
case 2: /* sparse */
get_first(rhs,"(");
str=get_next(",");
ntype=SPARSE;
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(",");
strcpy(indname,str);
iind=find_lookup(indname);
if(iind<0){
printf("in network %s, %s is not a table \n",
name,indname);
return 0;
}
str=get_next(")");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
my_net[ind].values=(double *)malloc((ntot+1)*sizeof(double));
init_net(my_net[ind].values,ntot);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].index=my_table[iind].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].n=ntot;
my_net[ind].ncon=ncon;
printf(" Added sparse %s len=%d x %d using %s var[%d] and %s\n",
name,ntot,ncon,wgtname,ivar,indname );
return 1;
break;
case 3: /* convolution */
get_first(rhs,"(");
str=get_next(",");
ntype=-1;
if(str[0]=='E')ntype=FCONVE;
if(str[0]=='0'||str[0]=='Z')ntype=FCONV0;
if(str[0]=='P')ntype=FCONVP;
if(ntype==-1){
printf(" No such convolution type %s \n",str);
return 0;
}
str=get_next(",");
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(",");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
str=get_next(",");
strcpy(root2name,str);
ivar2=get_var_index(root2name);
if(ivar2<0){
printf(" In %s , %s is not valid variable\n",
name,root2name);
return 0;
}
str=get_next(")");
strcpy(fname,str);
sprintf(junk,"%s(%s,%s)",fname,rootname,root2name);
if(add_expr(junk,my_net[ind].f,&elen)){
printf(" bad function %s \n",fname);
return 0;
}
my_net[ind].values=(double *)malloc((ntot+1)*sizeof(double));
init_net(my_net[ind].values,ntot);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].root2=ivar2;
my_net[ind].n=ntot;
my_net[ind].ncon=ncon;
printf(" Added net %s type %d len=%d x %d using %s %s(var[%d],var[%d]) \n",
name,ntype,ntot,ncon,wgtname,fname,ivar,ivar2);
return 1;
break;
case 4: /* sparse */
get_first(rhs,"(");
str=get_next(",");
ntype=FSPARSE;
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(",");
strcpy(indname,str);
iind=find_lookup(indname);
if(iind<0){
printf("in network %s, %s is not a table \n",
name,indname);
return 0;
}
str=get_next(",");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
str=get_next(",");
strcpy(root2name,str);
ivar2=get_var_index(root2name);
if(ivar2<0){
printf(" In %s , %s is not valid variable\n",
name,root2name);
return 0;
}
str=get_next(")");
strcpy(fname,str);
sprintf(junk,"%s(%s,%s)",fname,rootname,root2name);
if(add_expr(junk,my_net[ind].f,&elen)){
printf(" bad function %s \n",fname);
return 0;
}
my_net[ind].values=(double *)malloc((ntot+1)*sizeof(double));
init_net(my_net[ind].values,ntot);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].index=my_table[iind].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].root2=ivar2;
my_net[ind].n=ntot;
my_net[ind].ncon=ncon;
printf(" Sparse %s len=%d x %d using %s %s(var[%d],var[%d]) and %s\n",
name,ntot,ncon,wgtname,fname,ivar,ivar2,indname );
return 1;
break;
case 5: /* fft convolution */
get_first(rhs,"(");
str=get_next(",");
ntype=-1;
/* if(str[0]=='E')ntype=CONVE; */
if(str[0]=='0'||str[0]=='Z')ntype=FFTCON0;
if(str[0]=='P')ntype=FFTCONP;
if(ntype==-1){
printf(" No such fft convolution type %s \n",str);
return 0;
}
str=get_next(",");
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
ntab=get_lookup_len(iwgt);
if(type==FFTCONP&&ntab<ntot){
printf(" In %s, weight is length %d < %d \n",name,ntab,ntot);
return 0;
}
if(type==FFTCON0&&ntab<(2*ntot)){
printf(" In %s, weight is length %d < %d \n",name,ntab,2*ntot);
return 0;
}
str=get_next(")");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
if(ntype==FFTCON0)
ncon=2*ntot;
else
ncon=ntot;
my_net[ind].fftr=(double *)malloc((ncon+2)*sizeof(double));
my_net[ind].ffti=(double *)malloc((ncon+2)*sizeof(double));
my_net[ind].dr=(double *)malloc((ncon+2)*sizeof(double));
my_net[ind].di=(double *)malloc((ncon+2)*sizeof(double));
my_net[ind].iwgt=iwgt;
my_net[ind].values=(double *)malloc((ntot+1)*sizeof(double));
init_net(my_net[ind].values,ntot);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].n=ntot;
my_net[ind].ncon=ncon;
update_fft(ind);
printf(" Added net %s type %d len=%d x %d using %s var[%d] \n",
name,ntype,ntot,ncon,wgtname,ivar);
return 1;
break;
case 6: /* MMULT ntot=n,ncon=m */
get_first(rhs,"(");
str=get_next(",");
ntype=MMULT;
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(")");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
my_net[ind].values=(double *)malloc((ncon+1)*sizeof(double));
init_net(my_net[ind].values,ncon);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].n=ncon;
my_net[ind].ncon=ntot;
printf(" Added mmult %s len=%d x %d using %s var[%d]\n",
name,ntot,ncon,wgtname,ivar,indname );
return 1;
break;
case 7: /* FMMULT */
get_first(rhs,"(");
str=get_next(",");
ntype=FMMULT;
ntot=atoi(str);
if(ntot<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
ncon=atoi(str);
if(ncon<=0){
printf(" %s must be positive int \n",str);
return 0;
}
str=get_next(",");
strcpy(wgtname,str);
iwgt=find_lookup(wgtname);
if(iwgt<0){
printf("in network %s, %s is not a table \n",
name,wgtname);
return 0;
}
str=get_next(",");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
str=get_next(",");
strcpy(root2name,str);
ivar2=get_var_index(root2name);
if(ivar2<0){
printf(" In %s , %s is not valid variable\n",
name,root2name);
return 0;
}
str=get_next(")");
strcpy(fname,str);
sprintf(junk,"%s(%s,%s)",fname,rootname,root2name);
if(add_expr(junk,my_net[ind].f,&elen)){
printf(" bad function %s \n",fname);
return 0;
}
my_net[ind].values=(double *)malloc((ncon+1)*sizeof(double));
init_net(my_net[ind].values,ncon);
my_net[ind].weight=my_table[iwgt].y;
my_net[ind].type=ntype;
my_net[ind].root=ivar;
my_net[ind].root2=ivar2;
my_net[ind].n=ncon;
my_net[ind].ncon=ntot;
printf(" Added fmmult %s len=%d x %d using %s %s(var[%d],var[%d])\n",
name,ntot,ncon,wgtname,fname,ivar,ivar2);
return 1;
case 30:
/* interpolation array
z=INTERP(meth,n,root)
*/
get_first(rhs,"(");
str=get_next(",");
ivar=atoi(str);
my_net[ind].type=INTERP;
my_net[ind].iwgt=ivar;
str=get_next(",");
ivar=atoi(str);
if(ivar<1){
printf("Need more than 1 entry for interpolate\n");
return 0;
}
my_net[ind].n=ivar; /* # entries in array */
str=get_next(")");
strcpy(rootname,str);
ivar=get_var_index(rootname);
if(ivar<0){
printf(" In %s , %s is not valid variable\n",
name,rootname);
return 0;
}
my_net[ind].root=ivar;
printf("Added interpolator %s length %d on %s \n",name,my_net[ind].n,rootname);
return 1;
case 10:
/*
z=GILL(meth,rxn list)
e.g
z=GILL(meth,r{1-15})
GILL is different -
iwgt=evaluation method - 0 is standard
1 - tau-leap
root=number of reactions
values[0]=time of next reaction
values[1..root]=number of times this rxn took place
gcom contains list of all the fixed holding the reactions
*/
get_first(rhs,"(");
str=get_next(",");
ivar=atoi(str);
str=get_next(")");
my_net[ind].type=GILLTYPE;
my_net[ind].iwgt=ivar;
my_net[ind].gcom=(int *)malloc(1000*sizeof(int));
if(gilparse(str,my_net[ind].gcom,&ivar2)==0)
return 0;
my_net[ind].root=ivar2;
my_net[ind].n=ivar2+1;
my_net[ind].ncon=-1;
my_net[ind].values=(double *)malloc((ivar2+2)*sizeof(double));
printf("Added gillespie chain with %d reactions \n",ivar2);
return 1;
/* case 8:
get_first(rhs,"(");
str=get_next(",");
ntot=atoi(str);
str=get_next("{");
i=0;
elen=strlen(str);
while(1){
cc=str[i];
if(cc=='}'){junk[i]=0;
break;
}
junk[i]=cc;
i++;
if(i==elen){
printf("Illegal syntax for GROUP %s \n",str);
return 0;
}
}
printf("total=%d str=%s\n",ntot,junk);
return 0; */
}
return 0;
}
add_special_name(name,rhs)
char *name;
char *rhs;
{
if(is_network(rhs)){
printf(" netrhs = |%s| \n",rhs);
if(n_network>=MAXNET){
return;
}
strcpy(my_net[n_network].name,name);
add_net_name(n_network,name);
n_network++;
}
else
printf(" No such special type ...\n");
}
is_network(s)
char *s;
{
int i,n;
de_space(s);
strupr(s);
n=strlen(s);
if(s[0]=='C' &&s[1]=='O' &&s[2]=='N' && s[3]=='V')return 1;
if(s[0]=='S' &&s[1]=='P' &&s[2]=='A' && s[3]=='R')return 2;
if(s[0]=='F'&&s[1]=='C' &&s[2]=='O' &&s[3]=='N' && s[4]=='V')return 3;
if(s[0]=='F' && s[1]=='S' &&s[2]=='P' &&s[3]=='A' && s[4]=='R')return 4;
if(s[0]=='F' && s[1]=='F' && s[2]=='T' && s[3]=='C' )return 5;
if(s[0]=='M' &&s[1]=='M' &&s[2]=='U' && s[3]=='L')return 6;
if(s[0]=='F'&& s[1]=='M' &&s[2]=='M' &&s[3]=='U' && s[4]=='L')return 7;
if(s[0]=='G'&& s[1]=='I' &&s[2]=='L' &&s[3]=='L')return 10;
if(s[0]='I' && s[1]=='N' &&s[2]=='T' && s[3]=='E' && s[4]=='R')return INTERP;
/* if(s[0]=='G'&& s[1]=='R' && s[2]=='O' && s[3]=='U')return 8; */
return 0;
}
eval_all_nets()
{
int i;
for(i=0;i<n_network;i++)
evaluate_network(i);
}
evaluate_network(ind)
int ind;
{
int i,j,k,ij;
double sum,z;
int n=my_net[ind].n,*f;
int ncon=my_net[ind].ncon;
double *w,*y,*cc,*values,*yp;
int twon=2*n,root2=my_net[ind].root2;
cc=my_net[ind].index;
w=my_net[ind].weight;
values=my_net[ind].values;
y=&variables[my_net[ind].root];
switch(my_net[ind].type){
case INTERP: /* do nothing! */
break;
case GILLTYPE:
if(my_net[ind].ncon==-1&&my_net[ind].iwgt>0){
my_net[ind].weight=(double *)malloc(my_net[ind].root*NODE*sizeof(double));
make_gill_nu(my_net[ind].weight,NODE,my_net[ind].root,my_net[ind].values);
my_net[ind].ncon=0;
}
one_gill_step(my_net[ind].iwgt,my_net[ind].root,my_net[ind].gcom,my_net[ind].values);
break;
case CONVE:
for(i=0;i<n;i++){
sum=0.0;
for(j=-ncon;j<=ncon;j++){
k=abs(i+j);
if(k<twon){
if(k>=n)k=abs(twon-2-k);
sum+=(w[j+ncon]*y[k]);
}
}
values[i]=sum;
}
break;
case CONV0:
for(i=0;i<n;i++){
sum=0.0;
for(j=-ncon;j<=ncon;j++){
k=i+j;
if(k<n&&k>=0)
sum+=(w[j+ncon]*y[k]);
}
values[i]=sum;
}
break;
case CONVP:
for(i=0;i<n;i++){
sum=0.0;
for(j=-ncon;j<=ncon;j++){
k=((twon+i+j)%n);
sum+=(w[j+ncon]*y[k]);
}
values[i]=sum;
}
break;
case FFTCONP:
fft_conv(0,n,values,y,my_net[ind].fftr,my_net[ind].ffti,my_net[ind].dr,my_net[ind].di);
break;
case FFTCON0:
fft_conv(1,n,values,y,my_net[ind].fftr,my_net[ind].ffti,my_net[ind].dr,my_net[ind].di);
break;
case MMULT:
for(j=0;j<n;j++){
sum=0.0;
for(i=0;i<ncon;i++){
ij=j*ncon+i;
sum+=(w[ij]*y[i]);
}
values[j]=sum;
}
break;
case SPARSE:
for(i=0;i<n;i++){
sum=0.0;
for(j=0;j<ncon;j++){
ij=i*ncon+j;
k=(int)cc[ij];
if(k>=0&&k<n)
sum+=(w[ij]*y[k]);
}
values[i]=sum;
}
break;
/* f stuff */
case FCONVE:
f=my_net[ind].f;
yp=&variables[root2];
for(i=0;i<n;i++){
sum=0.0;
f[3]=(int)(&yp[i]);
for(j=-ncon;j<=ncon;j++){
k=abs(i+j);
if(k<twon){
if(k>=n)k=abs(twon-2-k);
f[1]=(int)(&y[k]);
z=evaluate(f);
sum+=(w[j+ncon]*z);
}
}
values[i]=sum;
}
break;
case FCONV0:
f=my_net[ind].f;
yp=&variables[root2];
for(i=0;i<n;i++){
sum=0.0;
f[3]=(int)(&yp[i]);
for(j=-ncon;j<=ncon;j++){
k=i+j;
if(k<n&&k>=0){
f[1]=(int)(&y[k]);
z=evaluate(f);
sum+=(w[j+ncon]*z);
}
}
values[i]=sum;
}
break;
case FCONVP:
f=my_net[ind].f;
yp=&variables[root2];
for(i=0;i<n;i++){
f[3]=(int)(&yp[i]);
sum=0.0;
for(j=-ncon;j<=ncon;j++){
k=((twon+i+j)%n);
f[1]=(int)(&y[k]);
z=evaluate(f);
sum+=(w[j+ncon]*z);
}
values[i]=sum;
}
break;
case FSPARSE:
f=my_net[ind].f;
yp=&variables[root2];
for(i=0;i<n;i++){
f[3]=(int)(&yp[i]);
sum=0.0;
for(j=0;j<ncon;j++){
ij=i*ncon+j;
k=(int)cc[ij];
if(k>=0&&k<n){
f[1]=(int)(&y[k]);
z=evaluate(f);
sum+=(w[ij]*z);
}
}
values[i]=sum;
}
break;
case FMMULT:
f=my_net[ind].f;
yp=&variables[root2];
for(j=0;j<n;j++){
f[3]=(int)(&yp[j]);
sum=0.0;
for(i=0;i<ncon;i++){
ij=j*ncon+i;
f[1]=(int)(&y[i]);
z=evaluate(f);
sum+=(w[ij]*z);
}
values[j]=sum;
}
break;
}
}
update_all_ffts()
{
int i;
for(i=0;i<n_network;i++)
if(my_net[i].type==FFTCON0||my_net[i].type==FFTCONP)
update_fft(i);
}
/*
tabular weights are of size 2k+1
and go from -k ... k
for FFT's
they are reordered as follows
fftr[i]=wgt[i+k] i = 0.. k
fftr[i+k]=wgt[i] i=1 .. k-1
*/
update_fft(int ind)
{
int i;
int dims[2];
double *w=my_net[ind].weight;
double *fftr=my_net[ind].fftr;
double *ffti=my_net[ind].ffti;
int n,n2;
int type=my_net[ind].type;
if(type==FFTCONP){
n=my_net[ind].n;
n2=n/2;
for(i=0;i<n;i++)ffti[i]=0.0;
for(i=0;i<=n2;i++)
fftr[i]=w[i+n2];
for(i=0;i<n2;i++)
fftr[n2+i+1]=w[i];
dims[0]=n;
fftn(1,dims,fftr,ffti,1,1.);
/* printf("index=%d n=%d n2=%d \n",ind,n,n2);
for(i=0;i<n;i++)
printf("(%g , %g)\n",fftr[i],ffti[i]); */
}
if(type==FFTCON0){
n=2*my_net[ind].n;
n2=n/2;
for(i=0;i<n;i++)ffti[i]=0.0;
for(i=0;i<=n2;i++)
fftr[i]=w[i+n2];
for(i=1;i<n2;i++)
fftr[n2+i]=w[i];
dims[0]=n;
fftn(1,dims,fftr,ffti,1,1.);
}
/* for(i=0;i<10;i++)printf("fftr,i=%g %g %g %g\n",fftr[i],ffti[i],fftr[n-1-i],ffti[n-1-i]); */
}
fft_conv(int it,int n,double *values,double *yy,double *fftr,double *ffti,double *dr,double *di)
{
int i;
int dims[2];
double x,y,mid;
int n2=2*n;
switch(it){
case 0:
dims[0]=n;
for(i=0;i<n;i++){
di[i]=0.0;
dr[i]=yy[i];
}
fftn(1,dims,dr,di,1,-2.0);
for(i=0;i<n;i++){
x=dr[i]*fftr[i]-di[i]*ffti[i];
y=dr[i]*ffti[i]+di[i]*fftr[i];
dr[i]=x;
di[i]=y;
}
fftn(1,dims,dr,di,-1,-2.0);
for(i=0;i<n;i++)
values[i]=dr[i];
return;
case 1:
dims[0]=n2;
for(i=0;i<n2;i++){
di[i]=0.0;
if(i<n)
dr[i]=yy[i];
else
dr[i]=0.0;
}
fftn(1,dims,dr,di,1,-2.0);
for(i=0;i<n2;i++){
x=dr[i]*fftr[i]-di[i]*ffti[i];
y=dr[i]*ffti[i]+di[i]*fftr[i];
dr[i]=x;
di[i]=y;
}
fftn(1,dims,dr,di,-1,-2.0);
for(i=0;i<n;i++)
values[i]=dr[i];
return;
}
}
/* parsing stuff to get gillespie code quickly */
gilparse(char *s,int *ind,int *nn)
{
int i=0,n=strlen(s);
char piece[50],b[20],bn[25],c;
int i1,i2,jp=0,f;
int k=0,iv;
int id,m;
printf("s=|%s|",s);
while(1){
c=s[i];
if(c==','||i>(n-1)){
piece[jp]=0;
if(g_namelist(piece,b,&f,&i1,&i2)==0){
printf("Bad gillespie list %s\n",s);
return 0;
}
if(f==0)
{
printf("added %s\n",b);
iv=get_var_index(b);
if(iv<0){
printf("No such name %s\n",b);
return 0;
}
ind[k]=iv;
k++;
}
else
{
printf("added %s{%d-%d}\n",b,i1,i2);
m=i2-i1+1;
for(id=0;id<m;id++){
sprintf(bn,"%s%d",b,id+i1);
iv=get_var_index(bn);
if(iv<0){
printf("No such name %s\n",bn);
return 0;
}
ind[k]=iv;
k++;
}
}
if(i>(n-1)){
*nn=k;
return 1;
}
jp=0;
}
else
{
piece[jp]=c;
jp++;
}
i++;
}
*nn=k;
return 1;
}
/* plucks info out of xxx{aa-bb} or returns string */
int g_namelist(char *s,char *root,int *flag,int *i1,int*i2)
{
int i,n=strlen(s),ir=-1,j=0;
char c,num[20];
*flag=0;
for(i=0;i<n;i++)
if(s[i]=='{')ir=i;
if(ir<0){
strcpy(root,s);
return 1;
}
for(i=0;i<ir;i++)
root[i]=s[i];
root[ir]=0;
*flag=1;
j=0;
for(i=ir+1;i<n;i++){
c=s[i];
if(c=='-')break;
num[j]=c;
j++;
}
if(i==n){
printf("Illegal syntax %s\n",s);
return 0;
}
num[j]=0;
*i1=atoi(num);
ir=i+1;
j=0;
for(i=ir;i<n;i++){
c=s[i];
if(c=='}')break;
num[j]=c;
j++;
}
num[j]=0;
*i2=atoi(num);
return 1;
}
syntax highlighted by Code2HTML, v. 0.9.1