gwr.mixed <- function(formula, data, regression.points, fixed.vars,intercept.fixed=FALSE, bw, diagnostic=T,
             kernel="bisquare", adaptive=FALSE, p=2, theta=0, longlat=F,dMat)
{
   ##Record the start time
  timings <- list()
  timings[["start"]] <- Sys.time()
  ###################################macth the variables
  this.call <- match.call()
  p4s <- as.character(NA)
  #####Check the given data frame and regression points
  #####Regression points
  if (missing(regression.points))
  {
  	rp.given <- FALSE
    regression.points <- data
    rp.locat<-coordinates(data)
    hatmatrix<-T
  }
  else
  {
    rp.given <- TRUE
    hatmatrix<-F
    if (is(regression.points, "Spatial"))
    {
       rp.locat<-coordinates(regression.points)
    }
    else if (is.numeric(regression.points) && dim(regression.points)[2] == 2)
       rp.locat<-regression.points
    else
      {
        warning("Output loactions are not packed in a Spatial object,and it has to be a two-column numeric vector")
        rp.locat<-dp.locat
      }
  }
  ##Data points{
  if (is(data, "Spatial"))
  {
    p4s <- proj4string(data)
    dp.locat<-coordinates(data)
    data <- as(data, "data.frame")
  }
  else
  {
    stop("Given regression data must be Spatial*DataFrame")
  }
    #########Distance matrix is given or not
  dp.n <- nrow(dp.locat)
  rp.n <- nrow(rp.locat)
  if (missing(dMat))
  {
      DM.given<-F
      DM1.given<-F
      if(dp.n + rp.n <= 10000)
      {
        dMat <- gw.dist(dp.locat=dp.locat, rp.locat=rp.locat, p=p, theta=theta, longlat=longlat)
        DM.given<-T
      }
  }
  else
  {
    DM.given<-T
    DM1.given<-T
    dim.dMat<-dim(dMat)
    if (dim.dMat[1]!=dp.n||dim.dMat[2]!=rp.n)
       stop("Dimensions of dMat are not correct")
  }
  ####################
  ######Extract the data frame
   mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data"), names(mf), 0L)

    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    mt <- attr(mf, "terms")
    y <- model.extract(mf, "response")
    x <- model.matrix(mt, mf)
    idx1 <- match("(Intercept)", colnames(x))
    if(!is.na(idx1))
      colnames(x)[idx1]<-"Intercept"
  #colnames(x)[1]<-"Intercept"
  if (missing(fixed.vars))
  {
    warning("No independent variables in the formula is specified as fixed terms!")
    if(!intercept.fixed)
       stop("Please use basic GWR function to calibrate this model")
  }
  else
  {
     if(intercept.fixed)
        fixed.vars <- c("Intercept", fixed.vars)
  }
  idx.fixed <- match(fixed.vars, colnames(x))
  x1 <- x[, -idx.fixed]
  x2<- x[, idx.fixed]
  if (!is.null(x1)) x1 <- as.matrix(x1, nrow = dp.n)
  if (!is.null(x2)) x2 <- as.matrix(x2, nrow = dp.n)
  colnames(x1) <- colnames(x)[-idx.fixed]
  colnames(x2) <- colnames(x)[idx.fixed]
  y <- as.matrix(y, nrow = dp.n)
  model <- gwr.mixed.2(x1, x2, y, dp.locat, out.loc=rp.locat, adaptive=adaptive, bw=bw,
                      kernel=kernel, p=p, theta=theta, longlat=longlat,dMat)
  res <- list()
   res$local <- model$local
   res$global <- apply(model$global,2,mean,na.rm=T)
   mgwr.df <- data.frame(model$local, model$global)
   colnames(mgwr.df) <- c(paste(colnames(x1), "L", sep="_"), paste(colnames(x2), "F", sep="_"))
   rownames(rp.locat)<-rownames(mgwr.df)
  griddedObj <- F
    if (is(regression.points, "Spatial"))
    {
        if (is(regression.points, "SpatialPolygonsDataFrame"))
        {
           polygons<-polygons(regression.points)
           #SpatialPolygons(regression.points)
           #rownames(mgwr.df) <- sapply(slot(polygons, "polygons"),
                              #  function(i) slot(i, "ID"))
           SDF <-SpatialPolygonsDataFrame(Sr=polygons, data=mgwr.df, match.ID=F)
        }
        else
        {
           griddedObj <- gridded(regression.points)
           SDF <- SpatialPointsDataFrame(coords=rp.locat, data=mgwr.df, proj4string=CRS(p4s), match.ID=F)
           gridded(SDF) <- griddedObj
        }
    }
    else
        SDF <- SpatialPointsDataFrame(coords=rp.locat, data=mgwr.df, proj4string=CRS(p4s), match.ID=F)
 #
#   if (is(regression.points, "SpatialPolygonsDataFrame"))
#    {
#       polygons<-polygons(regression.points)
#       #SpatialPolygons(regression.points)
#       rownames(mgwr.df) <- sapply(slot(polygons, "polygons"),
#                          function(i) slot(i, "ID"))
#       SDF <-SpatialPolygonsDataFrame(Sr=polygons, data=mgwr.df)
#    }
#    else
#       SDF <- SpatialPointsDataFrame(coords=rp.locat, data=mgwr.df, proj4string=CRS(p4s), match.ID=F)

   res$SDF <- SDF
   if (diagnostic)
   {
      gwr.fitted <- function(x,b) apply(x*b,1,sum)
      edf <- gwr.mixed.trace(x1, x2, y, dp.locat, adaptive=adaptive, bw=bw,
               kernel=kernel, p=p, theta=theta, longlat=longlat,dMat=dMat)
      model2 <- gwr.mixed.2(x1, x2, y, dp.locat, adaptive=adaptive, bw=bw,
                kernel=kernel, p=p, theta=theta, longlat=longlat,dMat=dMat)
      rss <- sum((y - gwr.fitted(model2$global, x2) - gwr.fitted(model2$local,x1))^2)
      n1 <- length(y)
      sigma.aic <- rss / n1
      aic <- log(sigma.aic*2*pi) + 1 + 2*(edf + 1)/(n1 - edf - 2)
      aic <- n1*aic
      res$aic <- aic
      res$df.used <- edf
      res$rss <- rss
   }
   GW.arguments<-list(formula=formula,rp.given=rp.given,hatmatrix=hatmatrix,bw=bw,
                       kernel=kernel,adaptive=adaptive, p=p, theta=theta, longlat=longlat,DM.given=DM1.given,diagnostic=diagnostic)
   res$GW.arguments <- GW.arguments
   res$this.call <- this.call
   timings[["stop"]] <- Sys.time()
   res$timings <- timings
   class(res) <- "mgwr"
   res}