# Coupled Damped Linear Oscillator using Latent Differnetial Equation Modeling
#
# Based on an article by Boker, Neale & Rausch (2004)
# Adapted by Pascal R. Deboeck from Mx script writen by Steven M. Boker
# 
# Last Revision: August 22, 2009
#_____________________________________________________________________________________________________

#Useful Bits
	rm(list=ls())
	library(OpenMx)
	
	Embed <- function(x,E,tau) {
		len <- length(x)
		out <- x[1:(len-(E*tau)+tau)]
		for(i in 2:E) { out <- cbind(out,x[(1+((i-1)*tau)):(len-(E*tau)+(i*tau))]) }
		return(out)
		}

#Provide the Data
	series1 <- sin(c(1:200)*.2)+rnorm(200,0,.3)
	series2 <- sin(c(1:200)*.1)+rnorm(200,0,.3)
	
	varNames <- c("a1","a2","a3","a4","a5","b1","b2","b3","b4","b5")
	inputData <- cbind(Embed(series1,5,3),Embed(series2,5,3))
	dimnames(inputData) <- list(NULL, varNames)
	
#Constants
	Number.occasions <- 5
	Number.indicators <- 1
	Number.subjects <- 121
	Number.time <- 2
	Number.variables <- Number.occasions*Number.indicators*2

#Mutable Matrices
	#A matrix
	A.info <- array(NA,dim=c(8,8,5))
	A.info[,,1] <- "F"
	A.info[,,2] <- 0
				##Free, Start, Name, Lower, Upper
	A.info[7,1,] <- c("T",-0.1,"eta1",NA,NA)
	A.info[7,2,] <- c("T",0.0,"zeta1",NA,NA)
	A.info[3,8,] <- c("T",0.0,"gamma1",NA,NA)
	A.info[8,4,] <- c("T",-0.1,"eta2",NA,NA)
	A.info[8,5,] <- c("T",0.0,"zeta2",NA,NA)
	A.info[6,7,] <- c("T",0.0,"gamma2",NA,NA)
	A.info[3,7,] <- c("F",1.0,"fixed",NA,NA)
	A.info[6,8,] <- c("F",1.0,"fixed",NA,NA)
	
	#S matrix
	S.info <- array(NA,dim=c(8,8,5))
	S.info[,,1] <- "F"
	S.info[,,2] <- 0
				##Free, Start, Name, Lower, Upper
	S.info[1,1,] <- c("T",0.01,"var x",0,NA)
	S.info[2,2,] <- c("T",0.01,"var dx",0,NA)
	S.info[3,3,] <- c("T",0.01,"var d2x",0,NA)
	S.info[4,4,] <- c("T",0.01,"var y",0,NA)
	S.info[5,5,] <- c("T",0.01,"var dy",0,NA)
	S.info[6,6,] <- c("T",0.01,"var d2y",0,NA)
	S.info[2,1,] <- c("T",0.0,"cov x-dx",NA,NA)
	S.info[1,2,] <- S.info[2,1,]	
	S.info[5,4,] <- c("T",0.0,"cov y-dy",NA,NA)
	S.info[4,5,] <- S.info[5,4,]
		
#Mx Matrices
	mA <- mxMatrix(	type="Full", free=A.info[,,1], labels=A.info[,,3], values=as.numeric(A.info[,,2]),
				lbound=as.numeric(A.info[,,4]),ubound=as.numeric(A.info[,,5]), name="mA")
	mS <- mxMatrix(	type="Symm", free=S.info[,,1], labels=S.info[,,3], values=as.numeric(S.info[,,2]),
				lbound=as.numeric(S.info[,,4]),ubound=as.numeric(S.info[,,5]), name="mS")

	time <- c(1:Number.occasions)-mean(c(1:Number.occasions))
	mK <- mxMatrix(	type="Full", values=cbind(time^0,time^1,(time^2)/2), name="mK")
	
	mI <- mxMatrix(	type="Iden", ncol=8, name="mI")
	mN <- mxMatrix(	type="Full", nrow=Number.occasions, ncol=3, name="mN")
	mQ <- mxMatrix(	type="Full", nrow=Number.occasions, ncol=2, name="mQ")
	mU <- mxMatrix(	type="Unit", nrow=Number.occasions, ncol=1, name="mU")
	mO <- mxMatrix(	type="Unit", nrow=1, ncol=1, name="mO")
	mP <- mxMatrix(	type="Full", free=c("F",rep("T",Number.indicators-1)), values=c(1,rep(0,Number.indicators-1)),
				nrow=Number.indicators, ncol=1, name="mP")
	mT <- mxMatrix(	type="Full", value=Number.time, nrow=1, ncol=1, name="mT")
	mE <- mxMatrix(	type="Diag", ncol=Number.variables, name="mE", free="T", values=.5)
	
	mB <- mxMatrix(	type="Full", free="T", nrow=1, ncol=Number.variables, name="mB", dimnames=list(NULL,varNames))
	#End Matrices;

#Mx Model
	CoupledModel <- mxModel("CoupledModel",
		mA, mS, mK, mI, mN, mQ, mU, mO, mP, mT, mE, mB,
		mxAlgebra(mP%x%(mK*(mU%x%cbind(mO,mT,mT*mT))),name="mM"),
		mxAlgebra(rbind(cbind(mM,mN,mQ),cbind(mN,mM,mQ)),name="mL"),
		mxAlgebra(solve(mI-mA)%&%mS,name="mR"),
		mxAlgebra(mL%&%mR+mE,name="expCovariance", dimnames=list(varNames,varNames)),
		mxData(observed=inputData,type="raw"),
		mxFIMLObjective(covariance="expCovariance",mean="mB")
		)

	
	EstModel <- mxRun(CoupledModel)