续上一篇：https://bbs.pinggu.org/thread-6038753-1-1.html

# 编写plotTree函数，使用igraph包的方法，通过遍历特征表达式，实习二叉树的绘制
getSubTree <- function(treeExp, dstID=NULL){
  library(stringr)
  # 提取子集数量
  gSize <- str_count(treeExp, "<")
  # 随机选取一个子集，并返回提取下标
  if(is.null(dstID))
    dstID <- round(runif(1, 1, gSize))
  loc_start = str_locate_all(treeExp,"<")[[1]][dstID,1]
  sumPair = 0
  for(i in loc_start:nchar(treeExp)){
    tchar <- substr(treeExp, start = i, stop = i)
    if(tchar=="<")sumPair=sumPair+1
    else if (tchar==">")sumPair=sumPair-1
    if(sumPair==0)break
  }
  nodeContext=substr(treeExp,start = loc_start+1, stop = i-1)
  isLeaf=is.na(str_locate(nodeContext,"<")[1,1])
  names(isLeaf)=NULL
  return(list(subStr=substr(treeExp, start = loc_start, stop = i),
              start=loc_start, end=i, nodeContext=nodeContext, isLeaf=isLeaf))
}



# 该函数用于以递归形式获取边集、点集
# treeExp:特征表达式
# id:初始节点号，通常为1
getEdgeVR <- function(treeExp,id=1){
  library(stringr)
  arrows=NULL
  idcount <<- idcount+1
  idList <<- c(idList,idcount)
  if (is.na(str_locate(treeExp,"<")[1,1])){
    txt=treeExp
    verNames<<-c(verNames,txt)
    isLeaf<<-c(isLeaf,1)
    return(id)
  }else{
    s0=str_locate(treeExp,",")[1,1]
    txt=substr(treeExp,start = 3,stop = s0-1)
    verNames<<-c(verNames,txt)
  }
  isLeaf<<-c(isLeaf,0)
  subt=getSubTree(treeExp,1)
  arrows <- c(arrows,paste(id, "->", getEdgeVR(subt$nodeContext,id+1)))
  # 若有两个子节点
  if(substr(treeExp,start = subt$end+1,stop=subt$end+1)==","){
    subt2=getSubTree(substr(treeExp,start = subt$end+2,
                            stop=nchar(treeExp)),1)
    arrows <- c(arrows,paste(id,"->",getEdgeVR(subt2$nodeContext,idcount)))
  }
  return(arrows)
}



# 该函数基于特征表达式treeExp,绘制二叉树
plotTree <- function(treeExp){
  library(igraph)
  idcount<<-0
  idList<<-NULL
  verNames<<-NULL
  isLeaf<<-NULL
  arws=getEdgeVR(treeExp,1)
  p_vertices=data.frame(idList,verNames,isLeaf)
  p_edges <- NULL
  for(obj in arws){
    tmp <- strsplit(obj," -> ")[[1]]
    tmpN <- length(tmp)
    p_edges <- rbind(p_edges,data.frame(from=tmp[1:(tmpN-1)],
                                        to=tmp[2:tmpN]))
  }

  p_edges=p_edges[complete.cases(p_edges),]
  p_edges=unique(p_edges)
  p_vertices.color=rep("Turquoise",nrow(p_vertices))
  p_vertices.color[p_vertices$isLeaf==1]="Orange"
  gg <- graph.data.frame(d=p_edges,directed=F,vertices=p_vertices)
  plot(gg,layout=layout.reingold.tilford,
       vertex.label=as.character(p_vertices$verNames),
       vertex.label.dist=0,vertex.color=p_vertices.color,
       vertex.label.color='Maroon',vertex.label.cex=1.2)
}
# 基于plotTree绘制特征表达式的二叉树
# plotTree(out)

# 产生初始种群
# 在遗传编程方法构建特征，需要产生初始种群，种群有N个个体定义，N表示种群规模
# 最佳个体表示特征的最佳组合，因此可用给定的m作为基因数量创建个体，每个基因由随机生成的特征表达式表示

# 产生k个(种群规模)个体函数genIndividuals，ksubs表示每个个体对应的固定基因数量或者长度
# 其中getAdjust函数计算个体适应度

genIndividuals <- function(k, ksubs, nMax=10)
  {
  individuals=NULL
  adjusts = NULL
  for (i in 1:k)
    {
    # 每个个体都从数据集的特征中产生表达树，并组合成个体
    singleTerms <- NULL
    for(j in 1:ksubs)
      {
      singleTerms <- c(singleTerms, randomGetTree("vdata",vfeatures))
    }
    individuals <- rbind(individuals,singleTerms)
    adjusts=c(adjusts, getAdjust(singleTerms))
  }
  rownames(individuals)=NULL
  individuals=data.frame(individuals,stringsAsFactors = F)
  individuals$adjusts=adjusts
  return(individuals)
}


# 计算适应度：回归问题，通常计算交叉验证的误差平方和降低量作为适应度，
# 分类问题可以根据精度的提高量或者信息增益作为适应度
# 适应度越大表示个体对环境适应能力越强，就越可能携带优秀基因，即特征越有效

# treeExpArray是一组基因或特征表达式，算法依次将特征表达式转换为值

getAdjust <- function(treeExpArray)
  {
  tempData=NULL
  # for(i in 1:length(treeExpArray))
  #   {
  #   treeExp <- treeExpArray

  for(treeExp in treeExpArray)
    {
    feature=eval(parse(text = gsub('>', '', gsub('<', '', treeExp))))
    if(is.na(sd(feature)) || is.nan(sd(feature)) || sd(feature)==0)
      {
      feature=rep(0,NROW(feature))
    }
    tempData <- cbind(tempData,feature)
  }

  colnames(tempData)=paste("X", 1:NROW(treeExpArray),sep="")
  tempData=data.frame(tempData)
  tempData$Y=vdata$Y
  newErr <- 0
  for(i in 1:13){
    trainData=tempData[setdiff(1:13,i),]
    testData=tempData[i,]
    newfit <- lm(Y~., data=trainData)
    testData$newPred <- predict(newfit,testData)
    newErr <- newErr+sum(abs(testData$Y - testData$newPred)^2)
  }
  # stdErr是全局值，表示基于原始属性得到的误差平方和
  interval=stdErr - newErr
  if(interval<0) return(0)
  return(interval)
}


未完待续。。。