"cta.15" <- function(n=3,t=5000, cues = NULL,stim = NULL, act=NULL, inhibit=NULL, consume = NULL,ten = NULL, type="both",fast=2 ) {
#simulation of the CTA reparamaterization of the dynamics of action
#revised 6/10/22 to include a stimulation matrix
compare <- FALSE
if(n > 4){ colours <- rainbow(n)} else {colours <- c("black","blue", "red", "green") }

step <- .05
ten.start <- ten
act.start <- act

tendencies.m <- matrix(NA,ncol=t,nrow=n)
actions.m <- matrix(NA,ncol=t,nrow=n)
if(is.null(stim))  stim <- diag(1,n)  #added 6/10/22
if(is.null(cues)) {cues <- 2^(n-1:n)}
if(is.null(inhibit)) {inhibit <-  matrix(1,ncol=n,nrow=n)
                     diag(inhibit) <- .05}
if(n>1) {colnames(inhibit) <- rownames(inhibit) <- paste("A",1:n,sep="")}
if(is.null(consume) ) {consume <- diag(.05,ncol=n,nrow=n) }
excite <- diag(step,n)

#first run for time= t  to find the maximum values to make nice plots as well as to get the summary stats
if (is.null(ten.start)) {ten <- rep(0,n)} else {ten <- ten.start}
if(is.null(act.start) ) {act <- cues} else {act <- act.start}
maxact <- minact <-  minten <-  maxten <- 0
counts <- rep(0,n)
transitions <- matrix(0,ncol=n,nrow=n)
frequency <- matrix(0,ncol=n,nrow=n)
actions <- tendencies <- rep(0,n)
acts <- tends <-  rep(0,n)

colnames(frequency) <- paste("T",1:n,sep="")
rownames(frequency) <- paste("F",1:n,sep="")
names(tendencies) <- paste("T",1:n,sep="")
names(actions) <- paste("A",1:n,sep="")
old.act <- which.max(act)
for (i in 1:t) {
	ten <- cues %*% stim %*%  excite  + ten - act %*% excite %*% consume    #but where are the sensitivities?
	act <- ten %*% excite + act - act %*% excite %*% inhibit
    
	act[act<0] <- 0
	tendencies <- tendencies + ten
	actions <- actions + act
	maxact <- max(maxact,act)
	minact <- min(minact,act)
	maxten <- max(maxten,ten)
	minten <- min(minten,ten)
	
	which.act <- which.max(act)
	counts[which.act] <- counts[which.act]+1
	acts[which.act] <- acts[which.act] + act[which.act]
	tends[which.act] <- tends[which.act] + ten[which.act]
	transitions[old.act,which.act] <- transitions[old.act,which.act] + 1
	if(old.act!=which.act) { frequency[old.act,which.act] <- frequency[old.act,which.act] + 1
	                         frequency[which.act,which.act] <- frequency[which.act,which.act] +1}
	old.act <- which.act
	tendencies.m[,i] <- ten
	actions.m[,i] <- act
}
 yl <- range(tendencies.m)
 plot(tendencies.m[1,],ylim=yl,xlab="time",ylab="Tendency",col="black",typ="l",main="Action tendencies over time",lwd=2)
 for (j in 2:n) {points(tendencies.m[j,],lty=j,col=colours[j],typ="l",lwd=2)}
 
  yl <- range(actions.m)
 plot(actions.m[1,],ylim=yl,xlab="time",ylab="Action Strength",col="black",typ="l",main="Actions over time",lwd=2)
 for (j in 2:n) {points(actions.m[j,],lty=j,col=colours[j],typ="l",lwd=2)}
 
if(FALSE){
#now do various types of plots, depending upon the type of plot desired

plots <- 1
action <- FALSE 
#state diagrams plot two tendencies agaist each other over time
if (type!="none") {if (type=="state") { 
	op <- par(mfrow=c(1,1))
	if (is.null(ten.start)) {ten <- rep(0,n)} else {ten <- ten.start}
	if(is.null(act.start) ) {act <- cues} else {act <- act.start}	
	plot(ten[1],ten[2],xlim=c(minten,maxten),ylim=c(minten,maxten),col="black", 
        main="State diagram",xlab="Tendency 1", ylab="Tendency 2") 
    
	for (i in 1:t) {
			ten <- cues %*%  excite  + ten - act %*% excite %*% consume
			act <- ten %*% excite + act - act %*% excite %*% inhibit
			act[act<0] <- 0
			if(!(i %% fast)) points(ten[1],ten[2],col="black",pch=20,cex=.2)
			} 
	}  else {
			 
#the basic default is to plot action tendencies and actions in a two up graph			 
if(type=="both") {if(compare) {op <- par(mfrow=c(2,2))} else {op <- par(mfrow=c(2,1))}
		plots <- 2 } else {op <- par(mfrow=c(1,1))}



if (type=="action") {action <- TRUE} else {if(type=="tend" ) action <- FALSE} 

for (k in 1:plots) {
	if (is.null(ten.start)) {ten <- rep(0,n)} else {ten <- ten.start}
	if(is.null(act.start) ) {act <- cues} else {act <- act.start}	
	if(action )   plot(rep(1,n),act,xlim=c(0,t),ylim=c(minact,maxact),xlab="time",ylab="action", main="Actions over time") else plot(rep(1,n),ten,xlim=c(0,t),ylim=c(minten,maxten),xlab="time",ylab="action tendency",main="Action Tendencies over time") 
		for (i in 1:t) {
			ten <- cues %*%  excite  + ten - act %*% excite %*% consume
			act <- ten %*% excite + act - act %*% excite %*% inhibit
			act[act<0] <- 0
			if(!(i %% fast) ) {if( action) points(rep(i,n),act,col=colours,cex=.2) else points(rep(i,n),ten,col=colours,cex=.2) }}
action <- TRUE}
} }

}   # end of if FALSE
acts <- acts/counts
tends <- tends/counts
cta.df <- data.frame(cues=cues,time=round(counts/t,2),frequency = rowSums(frequency),tendencies = round(t(tendencies/t)),actions = round(acts))
results <- list(cues=cues,cta=cta.df,inihibition=inhibit,time = counts/t,frequency=frequency,tendencies = tendencies/t,actions = actions/t)
class(results) <- c("psych","cta")
return(results)
}

"print_psych.cta" <- 
function(x,digits=2,all=NULL) {
if(all) {unclass(x)
   print(x)} else {
cat("\n Cues Tendency Action model\n")
cat("Cue strength = ",x$cues)
cat("\nTime spent =  ",round(x$time,digits))
cat("\nAverage Tendency = ",round(x$tendencies))
cat("\nAverage Action strength = ",round(x$actions))
}}