1"cta" <- function(n=3,t=5000, cues = NULL, act=NULL, inhibit=NULL,expect = NULL, consume = NULL,tendency = NULL,tstrength=NULL, type="both",fast=2 ,compare=FALSE,learn=TRUE,reward=NULL) {
2#simulation of the CTA reparamaterization of the dynamics of action
3#note, this is all for operation in a constant environment - what is needed is wrap the world around this.
4#That is, actions actually change cues
5#this should be a function to do all the work. the rest is just stats and graphics.
6#MODEL (dif equations from paper)
7
8##here is the basic model
9tf <- function(tendency,cues,step,expect,act,consume) {
10	tf <-  tendency + cues %*%  step %*% expect   - act %*% step %*% consume }   #tendencies at t-t
11
12
13af <- function(act,tendency,step,tstrength,inhibit) {
14		af  <- tendency %*% step %*% tstrength + act - act %*% step %*% inhibit} #actions at t-t
15#the learning function
16ef <- function(expect,act,step,consume,reward) {if(learn) {
17		which.act <- which.max(act) #counting which act has won
18
19		if(old.act!=which.act) {
20		diag(temp) <- act %*% reward
21      expect <- expect + temp
22	                          #learn but only if a new action
23	  expect <- expect*1/tr(expect)    #standardize expectancies
24	  old.act <- which.act}}
25
26ef <- expect
27}
28
29
30##
31temp <- matrix(0,n,n)
32
33
34if(n > 4){ colours <- rainbow(n)} else {colours <- c("blue", "red", "black","green") }
35
36stepsize <- .05
37tendency.start <- tendency
38act.start <- act
39expect.start <- expect
40
41if(is.null(cues)) {cues <- 2^(n-1:n)}   #default cue str - cue vector represents inate strength of stim (cake vs peanut)
42
43if(is.null(inhibit)) {inhibit <-  matrix(1,ncol=n,nrow=n)   #loss matrix .05 on diag 1s off diag
44                     diag(inhibit) <- .05}
45
46if(is.null(tstrength)) tstrength <- diag(1,n)
47#if(is.null(tstrength)) tstrength <- diag(c(1,.5,.25))
48
49if(n>1) {colnames(inhibit) <- rownames(inhibit) <- paste("A",1:n,sep="")}
50if(is.null(consume) ) {consume <- diag(.03,ncol=n,nrow=n) }
51step <- diag(stepsize,n) #this is the n CUES x k tendencyDENCIES matrix for Cue-tendency excitation weights
52if(is.null(expect)) expect <- diag(1,n)  # a matrix of expectancies that cues lead to outcomes
53#first run for time= t  to find the maximum values to make nice plots as well as to get the summary stats
54if (is.null(tendency.start)) {tendency <- rep(0,n)} else {tendency <- tendency.start} #default tendency = 0
55if(is.null(act.start) ) {act <- cues} else {act <- act.start} #default actions = initial cues
56
57if(is.null(reward)) {reward <- matrix(0,n,n)
58                     diag(reward) <- c(rep(0,n-1),.05) } else {temp1 <- reward
59                     reward <- matrix(0,n,n)
60                     diag(reward) <- temp1 }
61
62#set up counters
63maxact <- minact <-  mintendency <-  maxtendency <- 0
64counts <- rep(0,n)
65transitions <- matrix(0,ncol=n,nrow=n)
66frequency <- matrix(0,ncol=n,nrow=n)
67colnames(frequency) <- paste("T",1:n,sep="")
68rownames(frequency) <- paste("F",1:n,sep="")
69old.act <- which.max(act)
70
71
72#MODEL (dif equations from paper)
73for (i in 1:t) {
74
75
76	tendency <- tf(tendency,cues,step,expect,act,consume)
77	act <- af (act,tendency,step,tstrength,inhibit)
78	act[act<0] <- 0
79
80#add learning
81	  expect <- ef(expect,act,step,consume,reward)
82
83
84#END OF MODEL
85
86
87
88#STATS
89	#calc max/min act/tendency
90	maxact <- max(maxact,act)
91	minact <- min(minact,act)
92	maxtendency <- max(maxtendency,tendency)
93	mintendency <- min(mintendency,tendency)
94
95	#count
96	which.act <- which.max(act) #counting which act has won
97	counts[which.act] <- counts[which.act]+1 #time table
98	transitions[old.act,which.act] <- transitions[old.act,which.act] + 1 #frequency table
99	if(old.act!=which.act) { frequency[old.act,which.act] <- frequency[old.act,which.act] + 1
100	                         frequency[which.act,which.act] <- frequency[which.act,which.act] +1
101
102	                         }  #learn but only if a new action
103	old.act <- which.act
104}
105
106#PLOTS
107#now do various types of plots, depending upon the type of plot desired
108plots <- 1
109action <- FALSE
110#state diagrams plot two tendencydencies agaist each other over time
111if (type!="none") {if (type=="state") {
112	op <- par(mfrow=c(1,1))
113	if (is.null(tendency.start)) {tendency <- rep(0,n)} else {tendency <- tendency.start}
114	if(is.null(act.start) ) {act <- cues} else {act <- act.start}
115	plot(tendency[1],tendency[2],xlim=c(mintendency,maxtendency),ylim=c(mintendency,maxtendency),col="black",
116        main="State diagram",xlab="tendency 1", ylab="tendency 2")
117
118	for (i in 1:t) {
119			tendency <- tf(tendency,cues,step,expect,act,consume)
120	        act <- af (act,tendency,step,tstrength,inhibit)
121	        #expect <- ef(expect,act,step,consume)
122	act[act<0] <- 0
123			if(!(i %% fast)) points(tendency[1],tendency[2],col="black",pch=20,cex=.2)
124			}
125	}  else {
126
127#the basic default is to plot action tendencydencies and actions in a two up graph
128if(type=="both") {if(compare) {op <- par(mfrow=c(2,2))} else {op <- par(mfrow=c(2,1))}
129		plots <- 2 } else {op <- par(mfrow=c(1,1))}
130
131
132
133if (type=="action") {action <- TRUE} else {if(type=="tendencyd" ) action <- FALSE}
134
135for (k in 1:plots) {
136	if (is.null(tendency.start)) {tendency <- rep(0,n)} else {tendency <- tendency.start}
137	if(is.null(act.start) ) {act <- cues} else {act <- act.start}
138	if(is.null(expect.start)) {expect <- diag(1,n)} else {expect <- expect.start}   # a matrix of expectancies that cues lead to outcomes
139
140	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),tendency,xlim=c(0,t),ylim=c(mintendency,maxtendency),xlab="time",ylab="action tendency",main="Action tendencies over time")
141		for (i in 1:t) {
142		tendency <- tf(tendency,cues,step,expect,act,consume)
143	act <- af (act,tendency,step,tstrength,inhibit)
144
145	act[act<0] <- 0
146
147	###
148	maxact <- max(maxact,act)
149	minact <- min(minact,act)
150	maxtendency <- max(maxtendency,tendency)
151	mintendency <- min(mintendency,tendency)
152
153	#count
154	which.act <- which.max(act) #counting which act has won
155	counts[which.act] <- counts[which.act]+1 #time table
156	transitions[old.act,which.act] <- transitions[old.act,which.act] + 1 #frequency table
157	if(old.act!=which.act) { frequency[old.act,which.act] <- frequency[old.act,which.act] + 1
158	#frequency[which.act,which.act] <- frequency[which.act,which.act] +1
159	                       expect <- ef(expect,act,step,consume,reward)
160	                         }  #learn but only if a new action
161	old.act <- which.act
162
163	##
164
165
166
167			if(!(i %% fast) ) {if( action) points(rep(i,n),act,col=colours,cex=.2) else points(rep(i,n),tendency,col=colours,cex=.2) }}
168action <- TRUE}
169} }
170results <- list(cues=cues,expectancy=expect,strength=tstrength,inihibition=inhibit,consumation=consume,reinforcement=reward, time = counts,frequency=frequency, tendency=tendency, act=act)
171return(results)
172}