rm(list=ls())
require(geomorph)
nomall="Ai_SouthAf.TPS"
Lycan=readland.tps(nomall)
wid=10; hgt=4
quartz(width=wid, height=hgt)
old=par(mar=c(1,1,0.2,0.2))
gpa=gpagen(Lycan)
cat("plot gpa /n")
plot(gpa)
mtext("All",side=1,line=-2, adj=0.9, col='black',cex=2)
mals=readland.tps("Ai_SouthAf_mals.TPS")
nmal=length(mals[1,1,])
gpm=gpagen(mals)
quartz(width=wid, height=hgt)
old=par(mar=c(1,1,0.2,0.2))
cat("plot gpm /n")
plot(gpm)
mtext("Males",side=1,line=-2, adj=0.9, col='blue',cex=2)
fems=readland.tps("Ai_SouthAf_fems.TPS")
gpf=gpagen(fems)
nfem=length(fems[1,1,])
quartz(width=wid, height=hgt)
old=par(mar=c(1,1,0.2,0.2))
cat("plot gpf /n")
plot(gpf)
mtext("Females",side=1,line=-2, adj=0.9, col='red',cex=2)
fem=mshape(gpa$coords[,,1:nfem])
male=mshape(gpa$coords[,,(nfem+1):(nfem+nmal)])
sex=rep(NaN,nmal+nfem); sex[1:nfem]=1;sex[(nfem+1):(nmal+nfem)]=-1

quartz(width=wid, height=hgt)
old=par(mar=c(1,1,0.2,0.2))
cat("plot fem \n")
plot(fem, col='red',pch=20)
cat("add female landmark IDs \n")
text(fem,labels=1:length(fem[,1]) ,col='red',pos=2)
cat("add male points \n")
points(male, col='blue', pch=20)
outl=c(1:4,10,16,28,29,30,20,23,NaN,5:10,NaN,3,8,4,9,3,NA,18,11:16,NaN,17:20,NaN,18,19,22,18,NaN,11,13,12,19,NaN,1,5,17,21,22,24,25,26,24,NaN,2,7,15,27,28,29,14,NaN,13,30,NaN,22,23)
cat("add male lines \n")
lines(male[outl,], col='blue')
mtext(nomall,side=1,line=-2,adj=0.85)
mtext(paste(nmal,"males"),side=1,line=-1,adj=0.75,col='blue')
mtext(paste(nfem,"females"),side=1,line=-1,adj=0.95,col='red')
quartz()
boxplot(gpa$Csize~sex, main=paste(nomall,"Csize vs sex"))
mtext('male',side=1,line=3,col='blue', adj=0.25, cex=2)
mtext('female',side=1,line=3,col='red', adj=0.75, cex=2)
quartz()
hist(gpa$Csize)
hist(gpm$Csize, add=TRUE,col='blue')
hist(gpf$Csize, add=TRUE,col='red')
mtext(nomall)