/* Modified by Klaus Gebhardt, 1997 */

#include <stdio.h>
#include <math.h>
#include <X11/Xlib.h>
#include<X11/Xutil.h>

#include "xpplim.h"

#include "browse.h"

#define MAX_LEN_SBOX 25

double evaluate();
double drand48();


int (*rhs)();
extern float **storage;
extern int storind,FOUR_HERE;
extern int NODE,INFLAG,NEQ,NJMP,FIX_VAR,NMarkov;
extern char uvar_names[MAXODE][12];
float **my_adj;
int adj_len;
float **my_h;


extern char *adj_hint[];
extern Window info_pop;
extern char *info_message;
struct {
  int here,col0,ncol,colskip;
  int row0,nrow,rowskip; 
  float **data;
  char firstcol[11];
} my_trans;
   
int TRANPOSE_HERE=0;

double ADJ_EPS=1.e-8,ADJ_ERR=1.e-3;
int ADJ_MAXIT=6,ADJ_HERE=0,H_HERE=0,h_len,HODD_EV=0;
extern double DELTA_T,BOUND;
int *coup_fun[MAXODE];
char *coup_string[MAXODE];

extern int *my_ode[];
extern int NSYM,NSYM_START,NCON,NCON_START;
extern Window main_win;
extern int DCURY;

init_trans()
{
  my_trans.here=0;
  strcpy(my_trans.firstcol,uvar_names[0]);
  my_trans.ncol=2;
  my_trans.nrow=1;
  my_trans.rowskip=1;
  my_trans.colskip=1;
  my_trans.row0=1;
  my_trans.col0=2;
}

do_transpose()
{
 int i,status;
 char *n[]={"Column 1","NCols","ColSkip","Row 1","NRows","RowSkip"};
 char values[6][MAX_LEN_SBOX];
 sprintf(values[0],"%s",my_trans.firstcol);
 sprintf(values[1],"%d",my_trans.ncol);
 sprintf(values[2],"%d",my_trans.colskip);
 sprintf(values[3],"%d",my_trans.row0);
 sprintf(values[4],"%d",my_trans.nrow);
 sprintf(values[5],"%d",my_trans.rowskip);
 if(my_trans.here){
   
   for(i=0;i<=my_trans.nrow;i++)free(my_trans.data[i]);
   free(my_trans.data);
   my_trans.here=0;
   data_back();
 }
 status=do_string_box(6,6,1,"Transpose Data",n,values,33);
 if(status!=0){
   find_variable(values[0],&i);
   if(i>-1)
     my_trans.col0=i+1;
   else
     {
       err_msg("No such columns");
       return 0;
     }
   strcpy(my_trans.firstcol,values[0]);
   i=atoi(values[4]);
   if(i>=NEQ)i=NEQ-1;
   my_trans.nrow=i;
   my_trans.ncol=atoi(values[1]);
   my_trans.colskip=atoi(values[2]);
   my_trans.row0=atoi(values[3]);
   my_trans.rowskip=atoi(values[5]);
   return (create_transpose());
 }
 return 0; 
     
 
}
 
create_transpose()
{
  int i,j;
  int inrow,incol;
  my_trans.data=(float **)malloc(sizeof(float *)*(NEQ+1));
  for(i=0;i<=my_trans.nrow;i++)
    my_trans.data[i]=(float *)malloc(sizeof(float)*my_trans.ncol);
  for(i=my_trans.nrow+1;i<=NEQ;i++)my_trans.data[i]=storage[i];
  for(j=0;j<my_trans.ncol;j++)
    my_trans.data[0][j]=j+1;

  
  for(i=0;i<my_trans.ncol;i++){
    incol=my_trans.col0-1+i*my_trans.colskip;
    if(incol>NEQ)
      incol=NEQ;
    for(j=0;j<my_trans.nrow;j++){
      inrow=my_trans.row0+j*my_trans.rowskip;
      if(inrow>storind)
	inrow=storind;
      my_trans.data[j+1][i]=storage[incol][inrow];
    }
  }
  
   my_browser.data=my_trans.data;
   my_browser.col0=1;
   refresh_browser(my_trans.ncol);
   my_trans.here=1;
   return 1;
}
			 

alloc_h_stuff()
{
  int i;
 for(i=0;i<NODE ;i++){
   coup_fun[i]=(int *)malloc(100*sizeof(int));
   coup_string[i]=(char *)malloc(80);
   strcpy(coup_string[i],"0");
 }
}
 

data_back()
{
 FOUR_HERE=0;
 my_browser.data=storage;
 my_browser.col0=1;
 refresh_browser(storind);
}

adj_back()
{
 if(ADJ_HERE){
   my_browser.data=my_adj;

   my_browser.col0=1;
   refresh_browser(adj_len);
 }
}

h_back()
{
 if(H_HERE){
   my_browser.data=my_h;

   my_browser.col0=1;
   refresh_browser(h_len);
 }
}

make_adj()
{
 Window temp=main_win;
 static char *n[]={"(N)ew adj","(M)ake H","(A)djoint","(O)rbit","(H)fun"};
 static char key[]="nmaoh";
 char ch;
 int i;
 ch=(char)pop_up_list(&temp,"Adjoint",n,key,5,10,0,10,10*DCURY+8,adj_hint,
		      info_pop,info_message);
 switch(ch){
 case 'n': 
   new_adjoint();
   break;
 case 'm':
   new_h_fun();
   break;
 case 'a':
   adj_back();
   break;
 case 'o':
   data_back();
   break;
 case 'h':
   h_back();
   break;
 
 }
 }

  
 new_h_fun()
{
 int i,n=2;
 if(!ADJ_HERE){
   err_msg("Must compute adjoint first!");
   return;
 }
  if(storind!=adj_len){
     err_msg("incompatible data and adjoint");
     return;
   }
 if(H_HERE){
   free(my_h[0]);
   free(my_h[1]);
   if(HODD_EV){
     free(my_h[2]);
     free(my_h[3]);
   }
   free(my_h);
   H_HERE=0;
   HODD_EV=0;
 }
   if(NEQ>2){
     HODD_EV=1;
     n=4;
   }
   h_len=storind;
   data_back(); 
   my_h=(float **)malloc(sizeof(float*)*(NEQ+1));
   for(i=0;i<n;i++)my_h[i]=(float *)malloc(sizeof(float)*h_len);
   for(i=n;i<=NEQ;i++)my_h[i]=storage[i];
   if(make_h(storage,my_adj,my_h,h_len,DELTA_T*NJMP,NODE )){
     H_HERE=1;
     h_back();
   }
 ping();
  
}

make_h(orb,adj,h,nt,dt,node)
float **orb,**adj,**h;
double dt;
int node,nt;
{

 int i,j,rval=0;
 float sum;
 double z;
 int n0=node+1+FIX_VAR,k2,k;
 char name[30];
 for(i=0;i<NODE ;i++){
   sprintf(name,"Coupling for %s eqn:",uvar_names[i]);
   new_string(name,coup_string[i]);
   if(add_expr(coup_string[i],coup_fun[i],&j)){
     err_msg("Illegal formula");
     goto bye;
   }
 }
 /*  formulae are fine .. lets do it ... */
 for(j=0;j<nt;j++)  /* j is phi variable  */
   {
     sum=0.0;
     
     for(k=0;k<nt;k++){
         k2=k+j;
       if(k2>=nt)k2=k2-nt+1;
       for(i=0;i<node;i++){
	 set_ivar(i+1,(double)orb[i+1][k]);
         set_ivar(i+n0+1,(double)orb[i+1][k2]);
       }
       z=0.0;
       

       for(i=0;i<node;i++){
	
	 z=evaluate(coup_fun[i]);
	
	 sum=sum+(float)z*adj[i+1][k];
       }
	
     }
     my_h[0][j]=orb[0][j];
     my_h[1][j]=sum*dt;
   }
   if(HODD_EV){
     for(k=0;k<nt;k++){
       k2=nt-k-1;
       my_h[2][k]=.5*(my_h[1][k]-my_h[1][k2]);
       my_h[3][k]=.5*(my_h[1][k]+my_h[1][k2]);
     }
   }
   rval=1;
      
 bye:
  NSYM=NSYM_START;
  NCON=NCON_START;
  return(rval);

}
 
   
 
   

new_adjoint()
{
 int i,n=NODE +1;
 if(ADJ_HERE){
   data_back();
   for(i=0;i<n;i++)free(my_adj[i]);
   free(my_adj);
   ADJ_HERE=0;
 }
 adj_len=storind;
 my_adj=(float **)malloc((NEQ+1)*sizeof(float *));
 for(i=0;i<n;i++)my_adj[i]=(float *)malloc(sizeof(float)*adj_len);
 for(i=n;i<=NEQ;i++)my_adj[i]=storage[i];
 if(adjoint(storage,my_adj,adj_len,DELTA_T*NJMP,ADJ_EPS,ADJ_ERR,ADJ_MAXIT,NODE )){
   ADJ_HERE=1;;
 adj_back();
 }
 ping();
}
   







/*    ADJOINT ROUTINE
 *
 *
      This assumes that you have already computed the periodic orbit
	and have stored in in an array **orbit
      including time in the first column
     
     
      The righthand sides of the equations are
	rhs(t,y,yp,n)
      and the coupling function for ``H'' functions is
	couple(y,yhat,f,n)

	where yhat is presynaptic and y is postynaptic
     variable.  f returns the coupling vector.

    adjoint is the same size as orbit and when returned has 
    t in the first column.
  */

 adjoint(orbit,adjnt,nt,dt,eps,minerr,maxit,node)
 float **orbit,**adjnt;
 double dt,eps,minerr;
 int nt,node,maxit;
 {
  double **jac,*yold,ytemp,*fold,*fdev;
  double *yprime,*work;
  double t,prod,sup,del;
  int i,j,k,l,k2,rval=0,iout;
  int n2=node*node;
  double error;
   
   work = (double *)malloc((n2+4*node)*sizeof(double));
   yprime = (double *)malloc(node*sizeof(double));
   yold=(double *)malloc(node*sizeof(double));
   fold=(double *)malloc(node*sizeof(double));
   fdev=(double *)malloc(node*sizeof(double));
  jac = (double **)malloc(n2*sizeof(double *));
  
  for(i=0;i<n2;i++)
  {
   jac[i]=(double *)malloc(nt*sizeof(double));
       if(jac[i]==NULL){
       err_msg("Insufficient storage");
	return(0);
	}
   }
  
  /*  Now we compute the
	transpose time reversed jacobian  --  this is complex !! */
  for(k=0;k<nt;k++){
	l=nt-1-k;  /* reverse the limit cycle  */
	for(i=0;i<node;i++)yold[i]=(double)orbit[i+1][l];
        rhs(0.0,yold,fold,node);
	for(j=0;j<node;j++){
		ytemp=yold[j];
		del=eps*fabs(ytemp);
		if(del<eps)del=eps;
		
		yold[j]+=del;
		rhs(0.0,yold,fdev,node);
		yold[j]=ytemp;
		for(i=0;i<node;i++)
			jac[i+node*j][k]=(fdev[i]-fold[i])/del;
		
	      }
 
      }
 
  /* now we iterate to get a good adjoint using implicit Euler's method */
  ytemp=0.0;
 for(i=0;i<node;i++){
	yold[i]=1.+.01*(drand48()-.5); /* random initial data */
	
	ytemp+=fabs(yold[i]);
	}
  for(i=0;i<node;i++){
    yold[i]=yold[i]/ytemp;
    fdev[i]=yold[i];
  }
	
  printf("%f %f \n",yold[0],yold[1]);
      
 for(l=0;l<maxit;l++){
	for(k=0;k<nt-1;k++){
		k2=k+1;
		if(k2>=nt)k2=k2-nt;
		if(step_eul(jac,k,k2,yold,work,node,dt)==0){
		  rval=0;
		  goto bye;
		} 
		/* rk_interp(jac,k,k2,yold,work,node,dt,5); */
	      }
	ytemp=0.0;
	error=0.0;
	
         for(i=0;i<node;i++){
	  if(fabs(yold[i])>BOUND){
	    
	     rval=0;
	  err_msg("Out of bounds");
	  goto bye;
	  }
	error+=fabs(yold[i]-fdev[i]);
	ytemp+=fabs(yold[i]);
	
	}
	
	 
	
        for(i=0;i<node;i++){ yold[i]=yold[i]/ytemp;
			     fdev[i]=yold[i];
			   }
	printf("%f %f \n",yold[0],yold[1]);
        printf("err=%f \n",error);
	if(error<minerr)break; /*  exit if error small   */
 }
 /*  onelast time to compute the adjoint  */
 prod=0.0;   /* for normalization   */
 t=0.0;
 for(k=0;k<nt;k++){
        l=nt-k-1;
	t+=dt; 
        for(i=0;i<node;i++)fdev[i]=(double)orbit[i+1][l];
	rhs(0.0,fdev,yprime,node);
	for(j=0;j<node;j++){
	adjnt[j+1][l]=(float)yold[j];
	prod+=yold[j]*yprime[j]*dt;
      }
	k2=k+1;
	if(k2>=nt)k2-=nt;
	 if(step_eul(jac,k,k2,yold,work,node,dt)==0){
	  rval=0;
	  goto bye;
	    }	 
	/* rk_interp(jac,k,k2,yold,work,node,dt,5); */
	

      }

        	 
	prod=prod/t;
  for(k=0;k<nt;k++){
     for(j=0;j<node;j++)adjnt[j+1][k]=adjnt[j+1][k]/(float)prod;
     adjnt[0][k]=orbit[0][k];
   }
  rval=1;

 bye: 
   free(work);  
   free(yprime);
   free(yold);
   free(fold);
   free(fdev);
   for(i=0;i<n2;i++)
   free(jac[i]);
   free(jac); 
   return(rval);
 }
 
 

eval_rhs(jac,k1,k2,t,y,yp,node)
     double t,**jac,*y,*yp;
     int node,k1,k2;
{
  int i;
  int j;
  for(j=0;j<node;j++){
    yp[j]=0.0;
    for(i=0;i<node;i++)
 yp[j]=yp[j]+(jac[i+j*node][k1]*(1.0-t)+jac[i+j*node][k2]*t)*y[i];
  }
}

rk_interp(jac,k1,k2,y,work,neq,del,nstep)
double *y,del,*work,**jac;
int neq,k1,k2,nstep;
{
 int i,j;
 double *yval[3],dt=del/nstep;
 double t=0.0,t1,t2;
 yval[0]=work;
 yval[1]=work+neq;
 yval[2]=work+neq+neq;
 for(j=0;j<nstep;j++)
 {
  eval_rhs(jac,k1,k2,t/del,y,yval[1],neq);
  for(i=0;i<neq;i++)
  {
   yval[0][i]=y[i]+dt*yval[1][i]/6.00;
   yval[2][i]=y[i]+dt*yval[1][i]*0.5;
  }
  t1=t+.5*dt;
  eval_rhs(jac,k1,k2,t1/del,yval[2],yval[1],neq);
  for(i=0;i<neq;i++)
  {
   yval[0][i]=yval[0][i]+dt*yval[1][i]/3.00;
   yval[2][i]=y[i]+.5*dt*yval[1][i];
  }
  eval_rhs(jac,k1,k2,t1/del,yval[2],yval[1],neq);
  for(i=0;i<neq;i++)
  {
   yval[0][i]=yval[0][i]+dt*yval[1][i]/3.000;
   yval[2][i]=y[i]+dt*yval[1][i];
  }
  t2=t+dt;
  eval_rhs(jac,k1,k2,t2/del,yval[2],yval[1],neq);
  for(i=0;i<neq;i++)y[i]=yval[0][i]+dt*yval[1][i]/6.00;
  t=t2;
 }
return(1);  
}
 
	 
step_eul(jac,k,k2,yold,work,node,dt)
double *work,*yold,**jac,dt;
int k,k2,node;
{

int j,i,n2=node*node,info;
int ipvt[MAXODE];
double *mat,*fold;
fold=work;
mat=work+node;

  for(j=0;j<node;j++){
    fold[j]=0.0;
    for(i=0;i<node;i++)
      fold[j]=fold[j]+jac[i+j*node][k]*yold[i];
  }
  for(j=0;j<node;j++)yold[j]=yold[j]+.5*dt*fold[j];
  for(i=0;i<n2;i++)mat[i]=-jac[i][k2]*dt*.5;
  for(i=0;i<node;i++)mat[i+i*node]=1.+mat[i+i*node];
  sgefa(mat,node,node,ipvt,&info);
if(info!=-1){
  
  err_msg("Univertible Jacobian");
  return(0);
}
sgesl(mat,node,node,ipvt,yold,0);
return(1);
}
		
		











