R/cgam.R

Defines functions ddb.penal_dexp ddb.penal db.penal varest print.anova.cgam anova.cgam cgam.pvz cgam.pv getbin getvars makebin Ord fwt_a fgr_a fwt fgr ddfun dfun pfun polr_getedf predict_polr dev_fun print.summary.cgam.polr summary.cgam.polr cgam.polr.fit cgam.polr cpfun sameside intri plotpersp.trispl s.conc.conc s.conv.conv tri_getedf trispl.fit irls plot.shapeselect CharToShape ShapeToChar ConstrALL ConstrGA make_form in.or.out shapes best.fit ShapeSelect fitted.cgam.polr fitted.trispl fitted.wps fitted.cgam coef.cgam.polr coef.trispl coef.wps coef.cgam makeamat_wps makedelta_wps s.decr.decr s.decr.incr s.incr.decr s.incr.incr wps.fit wps_getedf plotpersp.wps plotpersp.cgam plotpersp pred_del predict.cgam print.summary.wps print.summary.cgam summary.wps summary.cgam incconcave incconvex concave convex mondecr monincr make_delta_add makedelta umbrella.fun umbrella tree.fun tree s s.decr.conc s.decr.conv s.incr.conc s.incr.conv s.conc s.conv s.decr s.incr decr.conc incr.conc decr.conv incr.conv conc conv decr incr CicFamily CicFamily cgam.fit bmat.fun amat.fun cgam

Documented in best.fit cgam conc conv decr decr.conc decr.conv incr incr.conc incr.conv in.or.out Ord plotpersp predict.cgam s s.conc s.conc.conc s.conv s.conv.conv s.decr s.decr.conc s.decr.conv s.decr.decr s.decr.incr shapes ShapeSelect s.incr s.incr.conc s.incr.conv s.incr.decr s.incr.incr tree umbrella

######
#cgam#
######
cgam <- function(formula, nsim = 0, family = gaussian, cpar = 1.2, data = NULL, weights = NULL, sc_x = FALSE, sc_y = FALSE, pnt = TRUE, pen = 0, var.est = NULL)
{
	cl <- match.call()
    if (is.character(family)) 
    	family <- get(family, mode = "function", envir = parent.frame())
  	if (is.function(family)) 
     	family <- family()
  	if (is.null(family$family)) 
    	stop("'family' not recognized!")
    if (family$family == "ordered") {
		rslt <- cgam.polr(formula, data = data, weights = weights, family = family, nsim = nsim, cpar = cpar)
		rslt$call <- cl
		return (rslt)
	} else {
		labels <- NULL
  		mf <- match.call(expand.dots = FALSE)
		m <- match(c("formula", "data"), names(mf), 0L)
  		mf <- mf[c(1L, m)]
  		mf[[1L]] <- as.name("model.frame")
  		mf <- eval(mf, parent.frame())
 		ynm <- names(mf)[1]
  		mt <- attr(mf, "terms")
  		y <- model.response(mf, "any")
  		if (family$family == "binomial") {
			if (class(y) == "factor") {
				y <- ifelse(y == levels(y)[1], 0, 1)
			}
  		}
#additive
  		shapes1_add <- NULL; shapes2_add <- NULL; shapes_add <- NULL
  		xmat_add <- NULL; xnms_add <- NULL
  		tr <- NULL; pl <- NULL; umb <- NULL
  		tree.delta <- NULL; umbrella.delta <- NULL 
  		tid1 <- NULL; tid2 <- NULL; tpos2 <- 0
  		uid1 <- NULL; uid2 <- NULL; upos2 <- 0
  		nums_add <- NULL; ks_add <- list(); sps_add <- NULL; xid_add <- 1
  		zmat <- NULL; zid <- NULL; zid0 <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_param <- NULL; is_fac <- NULL; vals <- NULL; st <- 1; ed <- 1
  		ztb <- list(); iztb <- 1
		iadd <- 0
		varlist_add <- NULL
#warp
		warp.delta <- NULL
  		nums_wp <- NULL; ks_wp <- list(); ks_wps <- list(); sps_wp <- NULL
  		xmat_wp <- NULL; x1_wp <- NULL; x2_wp <- NULL; xnms_wp <- NULL; nks0_wp <- NULL; ks0_wp <- NULL; dc <- NULL
  		x1_wps <- NULL; x2_wps <- NULL; dcss <- list()
  		#zmat <- NULL; zid <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_fac <- NULL; is_param <- NULL; vals <- NULL; st <- 1; ed <- 1
		iwps <- 0
		varlist_wps <- NULL
#tri
		tri.delta <- NULL
  		nums_tri <- NULL; ks_tri <- list(); ks_tris <- list(); nks_tris <- list(); sps_tri <- NULL
  		xmat_tri <- NULL; x1_tri <- NULL; x2_tri <- NULL; xnms_tri <- NULL; nks0_tri <- NULL; ks0_tri <- NULL; cvs <- NULL	
  		x1_tris	<- NULL; x2_tris <- NULL; cvss <- list()
  		#zmat <- NULL; zid <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_fac <- NULL; is_param <- NULL; vals <- NULL; st <- 1; ed <- 1
  		#print (dim(mf))
		itri <- 0
		#print (head(mf))
  		for (i in 2:ncol(mf)) {
#additive
#print (attributes(mf[,i])$class)
		if (!is.null(attributes(mf[,i])$categ)) {
    		if (is.numeric(attributes(mf[,i])$shape) & (attributes(mf[,i])$categ == "additive")) {
       			labels <- c(labels, "additive")
       			shapes1_add <- c(shapes1_add, attributes(mf[,i])$shape)
       			xmat_add <- cbind(xmat_add, mf[,i])
       			xnms_add <- c(xnms_add, attributes(mf[,i])$nm)
       			nums_add <- c(nums_add, attributes(mf[,i])$numknots)
       			sps_add <- c(sps_add, attributes(mf[,i])$space)
       			ks_add[[xid_add]] <- attributes(mf[,i])$knots
       			xid_add <- xid_add + 1
    		}
    		if (is.character(attributes(mf[,i])$shape) & (attributes(mf[,i])$categ == "additive")) {
    			labels <- c(labels, "additive")
       			shapes2_add <- c(shapes2_add, attributes(mf[,i])$shape)
       			if (attributes(mf[,i])$shape == "tree") {
					pl <- c(pl, attributes(mf[,i])$pl)
					treei <- tree.fun(mf[,i], attributes(mf[,i])$pl)
					tree.delta <- rbind(tree.delta, treei)
					tpos1 <- tpos2 + 1
					tpos2 <- tpos2 + nrow(treei)
					tid1 <- c(tid1, tpos1)
					tid2 <- c(tid2, tpos2)
					tr <- cbind(tr, mf[,i])
       			}
       			if (attributes(mf[,i])$shape == "umbrella") {
					umbi <- umbrella.fun(mf[,i])
					umbrella.delta <- rbind(umbrella.delta, umbi)
					upos1 <- upos2 + 1
					upos2 <- upos2 + nrow(umbi)
					uid1 <- c(uid1, upos1)
					uid2 <- c(uid2, upos2)
					umb <- cbind(umb, mf[,i])	
				}
    		}
#warp
			if (is.character(attributes(mf[, i])$shape) & (attributes(mf[,i])$categ == "warp")) {
				iwps <- iwps + 1
				labels <- c(labels, rep(paste("warp", iwps, sep = "_"), 2))
    			sps_wp <- attributes(mf[, i])$space
				dcs <- attributes(mf[, i])$decreasing
				dcss[[iwps]] <- dcs
				ks0_wp <- attributes(mf[, i])$knots
				nks0_wp <- attributes(mf[, i])$numknots
				x1_wp <- (mf[, i])[, 1]
				x1_wps <- cbind(x1_wps, x1_wp)
				x2_wp <- (mf[, i])[, 2]
				x2_wps <- cbind(x2_wps, x2_wp)
				xmat_wp <- cbind(xmat_wp, mf[, i])
				xnms_wp <- c(xnms_wp, attributes(mf[, i])$name)
				#print (xnms_wp)
				#print (dim(xmat_wp))
				ans_warp <- makedelta_wps(x1t = x1_wp, x2t = x2_wp, m1_0 = nks0_wp[1], m2_0 = nks0_wp[2], k1 = ks0_wp$k1, k2 = ks0_wp$k2, space = sps_wp, decreasing = dcs)
				warp.delta0 <- ans_warp$delta
				k1 <- ans_warp$k1
				k2 <- ans_warp$k2
				ks_wp[[1]] <- k1
				ks_wp[[2]] <- k2
				ks_wps[[iwps]] <- ks_wp
				if (iwps > 1) {
					warp.delta0 <- warp.delta0[,-1]
				}
				varlist_wps <- c(varlist_wps, 1:ncol(warp.delta0)*0 + iwps)
				warp.delta <- cbind(warp.delta, warp.delta0)
    			#print (dim(warp.delta))
    			#print (ks_wps)
    		}
#tri
			if (is.character(attributes(mf[, i])$shape) & (attributes(mf[,i])$categ == "tri")) {
				itri <- itri + 1
				labels <- c(labels, rep(paste("tri", itri, sep = "_"), 2))
				sps_tri <- attributes(mf[, i])$space
				cvs <- attributes(mf[, i])$cvs
				cvss[[itri]] <- cvs
				ks0_tri <- attributes(mf[, i])$knots
				nks0_tri <- attributes(mf[, i])$numknots
				ks_tris[[itri]] <- ks0_tri
				nks_tris[[itri]] <- nks0_tri
				x1_tri <- (mf[, i])[, 1]
				x1_tris <- cbind(x1_tris, x1_tri)
				x2_tri <- (mf[, i])[, 2]
				x2_tris <- cbind(x2_tris, x2_tri)
				#xmat_tri <- cbind(x1_tri, x2_tri)
				xmat_tri <- cbind(xmat_tri, mf[, i])
				xnms_tri <- c(xnms_tri, attributes(mf[, i])$name)
    		}
    		#print (class((mf[, i])))
		}
    		if (is.null(attributes(mf[,i])$shape)) {
				if (!is.null(names(mf)[i])) {
	  				znms <- c(znms, names(mf)[i])
				}
        		if (!is.matrix(mf[,i])) {
          			zid <- c(zid, i)
	  				is_param <- c(is_param, TRUE)
          			if (is.factor(mf[,i])) {
	    				is_fac <- c(is_fac, TRUE)
	    				ch_char <- suppressWarnings(is.na(as.numeric(levels(mf[, i]))))
            			if (any(ch_char)) {
	      					vals <- c(vals, unique(levels(mf[, i]))[-1])
            			} else {
	      					vals <- c(vals, as.numeric(levels(mf[, i]))[-1])
	    				}
            			nlvs <- length(attributes(mf[,i])$levels)
	    				ed <- st + nlvs - 2 
	    				zid1 <- c(zid1, st)
	    				zid2 <- c(zid2, ed)
	    				st <- st + nlvs - 1
	    				zmat0 <- model.matrix(~ mf[, i])[, -1, drop = FALSE]
	    				zmat <- cbind(zmat, zmat0)
	    				ztb[[iztb]] <- mf[,i]
	    				iztb <- iztb + 1
          			} else {
	    				is_fac <- c(is_fac, FALSE)
            			zmat <- cbind(zmat, mf[, i])
	    				ztb[[iztb]] <- mf[,i]
            			iztb <- iztb + 1
	    				ed <- st
            			zid1 <- c(zid1, st)
	    				zid2 <- c(zid2, ed) 
	    				st <- st + 1
	    				vals <- c(vals, "")
         			}
       			} else {
	  				is_param <- c(is_param, FALSE)
          			is_fac <- c(is_fac, FALSE)
	  				zmat0 <- mf[, i]
	  				mat_cols <- ncol(zmat0)
	  				mat_rm <- NULL
	  				#rm_num <- 0
	  				for (irm in 1:mat_cols) {
       	  				if (all(round(diff(zmat0[, irm]), 8) == 0)) {
               				mat_rm <- c(mat_rm, irm)
          				}
   	  				}
	  				if (!is.null(mat_rm)) {
	  					zmat0 <- zmat0[, -mat_rm, drop = FALSE]
						#rm_num <- rm_num + length(mat_rm)
	  				}
	  				zmat <- cbind(zmat, zmat0)
         			ztb[[iztb]] <- mf[,i]
          			iztb <- iztb + 1
	  				vals <- c(vals, 1)
	  				zid <- c(zid, i)
	  				nlvs <- ncol(zmat0) + 1
	  				ed <- st + nlvs - 2
	  				zid1 <- c(zid1, st)
	  				zid2 <- c(zid2, ed)
	  				st <- st + nlvs - 1
      			}
    		}
  		}
  		dimnames(zmat)[[2]] <- NULL
  		if (family$family == "binomial" | family$family == "poisson" | family$family == "Gamma") {
    		wt.iter = TRUE
  		} else {wt.iter = FALSE}
  	  	#attr(xmat, "shape") <- shapes1
  	  	#print (labels)
  		if (any(labels == "additive") | !is.null(zmat)) {
  			xmat0_add <- xmat_add; shapes0_add <- shapes1_add; nums0_add <- nums_add; ks0_add <- ks_add; sps0_add <- sps_add; xnms0_add <- xnms_add; idx_s <- NULL; idx <- NULL
  			if (any(shapes1_add == 17)) {
    			kshapes <- length(shapes1_add)
    			obs <- 1:kshapes
    			idx_s <- obs[which(shapes1_add == 17)]; idx <- obs[which(shapes1_add != 17)]
    			xmat0_add[ ,1:length(idx_s)] <- xmat_add[ ,idx_s]
    			shapes0_add[1:length(idx_s)] <- shapes1_add[idx_s]
    			nums0_add[1:length(idx_s)] <- nums_add[idx_s]
    			sps0_add[1:length(idx_s)] <- sps_add[idx_s]
    			ks0_add[1:length(idx_s)] <- ks_add[idx_s]
    			xnms0_add[1:length(idx_s)] <- xnms_add[idx_s]
    			if (length(idx) > 0) {
      				xmat0_add[ ,(1 + length(idx_s)):kshapes] <- xmat_add[ ,idx]
      				shapes0_add[(1 + length(idx_s)):kshapes] <- shapes1_add[idx]
      				nums0_add[(1 + length(idx_s)):kshapes] <- nums_add[idx]
      				sps0_add[(1 + length(idx_s)):kshapes] <- sps_add[idx]
      				ks0_add[(1 + length(idx_s)):kshapes] <- ks_add[idx]
      				xnms0_add[(1 +length(idx_s)):kshapes] <- xnms_add[idx]
    			}
    			#xmat <- xmat0; nums <- nums0; ks <- ks0; sps <- sps0; xnms <- xnms0
  			}
  			shapes_add <- c(shapes1_add, shapes2_add)
  		}
  		#print (labels)
  		xnms <- c(xnms_wp, xnms_add, xnms_tri)
  		xmat <- cbind(xmat_wp, xmat_add, xmat_tri)
  		colnames(xmat) <- xnms
  		#boolz <- is.null(labels) & !is.null(zmat)
  		if (all(labels == "additive")) {
  			if (is.null(shapes_add)) {
    			nsim <- 0
  			}
  			ans <- cgam.fit(y = y, xmat = xmat0_add, zmat = zmat, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add, nsim = nsim, family = family, cpar = cpar, wt.iter = wt.iter, umbrella.delta = umbrella.delta, tree.delta = tree.delta, weights = weights, zid = zid, zid1 = zid1, zid2 = zid2, sc_x = sc_x, sc_y = sc_y, idx_s = idx_s, idx = idx, var.est = var.est, data = data)
  			if (!is.null(uid1) & !is.null(uid2)) {
    			uid1 <- uid1 + ans$d0 + ans$capm
    			uid2 <- uid2 + ans$d0 + ans$capm
  			}
  			if (!is.null(tid1) & !is.null(tid2)) {
    			tid1 <- tid1 + ans$d0 + ans$capm + ans$capu
    			tid2 <- tid2 + ans$d0 + ans$capm + ans$capu 
  			}
#new:
  			knots_add <- ans$knots
  			numknots_add <- ans$numknots
#new:
  			if (length(knots_add) > 0) {
    			names(knots_add) <- xnms_add
  			}
  			rslt <- list(etahat = ans$etahat, muhat = ans$muhat, vcoefs = ans$vcoefs, xcoefs = ans$xcoefs, zcoefs = ans$zcoefs, ucoefs = ans$ucoefs, tcoefs = ans$tcoefs, coefs = ans$coefs, cic = ans$cic, d0 = ans$d0, edf0 = ans$edf0, etacomps = ans$etacomps, y = y, xmat_add = xmat_add, zmat = zmat, ztb = ztb, tr = tr, umb = umb, tree.delta = tree.delta, umbrella.delta = umbrella.delta, bigmat = ans$bigmat, shapes = shapes_add, shapesx = shapes1_add, shapesx0 = shapes0_add, prior.w = weights, wt = ans$wt, wt.iter = ans$wt.iter, family = family, SSE0 = ans$sse0, SSE1 = ans$sse1, pvals.beta = ans$pvals.beta, se.beta = ans$se.beta, null_df = ans$df.null, df = ans$df, resid_df_obs = ans$resid_df_obs, edf = ans$df_obs, null_deviance = ans$dev.null, deviance = ans$dev, tms = mt, capm = ans$capm, capms = ans$capms, capk = ans$capk, capt = ans$capt, capu = ans$capu, xid1 = ans$xid1, xid2 = ans$xid2, tid1 = tid1, tid2 = tid2, uid1 = uid1, uid2 = uid2, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, nsim = nsim, xnms_add = xnms_add, xnms = xnms, ynm = ynm, znms = znms, is_param = is_param, is_fac = is_fac, knots = knots_add, numknots = numknots_add, sps = sps_add, ms = ans$ms, cpar = ans$cpar, pl = pl, idx_s = idx_s, idx = idx, xmat0 = ans$xmat2, knots0 = ans$knots2, numknots0 = ans$numknots2, sps0 = ans$sps2, ms0 = ans$ms2, phihat = ans$phihat, pvs = ans$pvs, s.edf = ans$s.edf, bstats = ans$bstats, pvsz = ans$pvsz, z.edf = ans$z.edf, fstats = ans$fstats, vh = ans$vh, kts.var = ans$kts.var)
  			rslt$call <- cl
  			class(rslt) <- "cgam"
		} else if (any(grepl("warp", labels, fixed = TRUE)) & !any(grepl("tri", labels, fixed = TRUE))) {
			nprs <- sum(grepl("warp", labels, fixed = TRUE)) / 2
			delta_add <- NULL
			np_add <- 0
			pb <- 0
			if (any(labels == "additive")) {
				#delta_add <- NULL
				if (length(shapes0_add) > 0) {
					ans_add <- make_delta_add(xmat = xmat0_add, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add) 
					delta_add <- ans_add$bigmat
#np_add: smooth only, conv and conc					
					np_add <- ans_add$np
					pb <- nrow(delta_add)
					varlist_add <- ans_add$varlist
				}
			}
			delta_ut <- NULL
			if (!is.null(umbrella.delta)) {
				delta_add <- rbind(delta_add, umbrella.delta)	
				delta_ut <- rbind(delta_ut, umbrella.delta)
			} 
			if (!is.null(tree.delta)) {
				delta_add <- rbind(delta_add, tree.delta)
				delta_ut <- rbind(delta_ut, tree.delta)
			}
			ans <- wps.fit(x1t = x1_wps, x2t = x2_wps, y = y, zmat = zmat, xmat_add = xmat_add, delta_add = delta_add, delta_ut = delta_ut, varlist_add = varlist_add, shapes_add = shapes_add, np_add = np_add, w = weights, pen = pen, pnt = pnt, cpar = cpar, decrs = dcss, delta = warp.delta, kts = ks_wps, wt.iter = wt.iter, family = family, nsim = nsim, nprs = nprs, idx_s = idx_s, idx = idx)
  			rslt <- list(ks_wps = ans$kts, muhat = ans$muhat, SSE1 = ans$sse1, SSE0 = ans$sse0, edfc = ans$edf, edf0 = ans$edf0, delta = ans$delta, y = y, zmat = ans$zmat, zmat_0 = ans$zmat_0, xmat_wp = xmat_wp, xmat_add = xmat_add, xmat = xmat, shapes = shapes_add, coefs = ans$coefs, zcoefs = ans$zcoefs, pvals.beta = ans$pz, se.beta = ans$sez, gcv = ans$gcv, xnms_wp = xnms_wp, xnms_add = xnms_add, xnms = xnms, xid_add = xid_add, znms = znms, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, ynm = ynm, decrs = dcss, tms = mt, is_param = is_param, is_fac = is_fac, d0 = ans$d0, pb = pb, pen = ans$pen, cpar = ans$cpar, wt.iter = wt.iter, cic = ans$cic, family = family, nprs_wps = nprs, varlist_wps = varlist_wps, coef_wp = ans$coef_wp, varlist_add = varlist_add, np_add = np_add, coef_add = ans$coef_add, coef_ut = ans$coef_ut, etacomps = ans$etacomps)
  			class(rslt) <- "wps"
  			if (!is.null(delta_add)) {
  				class(rslt) <- c("wps", "cgam")
  				if (min(which(labels == "additive")) < min(which(grepl("warp", labels, fixed = TRUE)))) {
  					class(rslt) <- c("cgam", "wps")
  				}
  			} 
  		} else if (any(grepl("tri", labels, fixed = TRUE))) {
  			nprs <- sum(grepl("tri", labels, fixed = TRUE)) / 2
  			nprs_wps <- sum(grepl("warp", labels, fixed = TRUE)) / 2
  			amat_wp <- NULL
  			dmat_wp <- NULL
			delta_add <- NULL
  			np_add <- 0
			pb <- 0
			if (any(labels == "additive")) {
				if (length(shapes0_add) > 0) {
					ans_add <- make_delta_add(xmat = xmat0_add, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add) 
					delta_add <- ans_add$bigmat
					np_add <- ans_add$np
					pb <- nrow(delta_add)
					varlist_add <- ans_add$varlist
				}
			}
			if (!is.null(umbrella.delta)) {
				delta_add <- rbind(delta_add, umbrella.delta)	
			} 
			if (!is.null(tree.delta)) {
				delta_add <- rbind(delta_add, tree.delta)
			}
			if (any(grepl("warp", labels, fixed = TRUE))) {
				ans_wps <- makeamat_wps(ks_wps, nprs_wps)
				amat_wp <- ans_wps$amat
				dmat_wp <- ans_wps$dmat
				valist_wps <- ans_wps$valist_wps
			}
  			ans <- trispl.fit(x1t = x1_tris, x2t = x2_tris, y = y, zmat = zmat, xmat_add = xmat_add, delta_add = delta_add, varlist_add = varlist_add, shapes_add = shapes_add, np_add = np_add, xmat_wp = xmat_wp, delta_wp = warp.delta, varlist_wps = varlist_wps, amat_wp = amat_wp, dmat_wp = dmat_wp, w = weights, lambda = pen, pnt = TRUE, cpar = cpar, cvss = cvss, delta = tri.delta, kts = ks_tris, nkts = nks_tris, wt.iter = wt.iter, family = family, nsim = nsim, nprs = nprs)
  			rslt <- list(muhat = ans$muhat, etahat = ans$etahat, trimat = ans$trimat, y = y, zmat = ans$zmat, capk = ans$capk, xmat_tri = xmat_tri, xmat_add = xmat_add, xmat_wp = xmat_wp, xmat = xmat, shapes = shapes_add, coefs = ans$thhat, zcoefs = ans$zcoefs, se.beta = ans$se.beta, pvals.beta = ans$pvals.beta, gcv = ans$gcv, edf = ans$edf, edf0 = ans$edf0, cic = ans$cic, sse0 = ans$sse0, sse1 = ans$sse1, family = family, nprs = nprs, nprs_wps = nprs_wps, varlist_wps = varlist_wps, varlist_add = varlist_add, varlist_tri = ans$varlist, xnms_tri = xnms_tri, xnms_add = xnms_add, xnms_wp = xnms_wp, xnms = xnms, znms = znms, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, ynm = ynm, decrs = dcss, cvss = cvss, tms = mt, is_param = is_param, is_fac = is_fac, d0 = ans$d0, pen = ans$pen, cpar = cpar, wt.iter = wt.iter, np_add = np_add, pb = pb, coef_add = ans$coef_add, coef_tri = ans$coef_tri, coef_wp = ans$coef_wp, etacomps = ans$etacomps, ks_wps = ks_wps, knots_lst = ans$knots_lst, trimat_lst = ans$trimat_lst, bmat_lst = ans$bmat_lst, capk_lst = ans$capk_lst)  				
  			class(rslt) <- "trispl"
  			if (!is.null(delta_add) & is.null(warp.delta)) {
  				if (labels[1] == "additive") {
  					class(rslt) <- c("cgam", "trispl")
  				} else {class(rslt) <- c("trispl", "cgam")} 				
  			}
  			if (!is.null(warp.delta) & is.null(delta_add)) {
  				#if (labels[1] == "warp") {
  				if (grepl("warp", labels[1], fixed = TRUE)) {
  					class(rslt) <- c("wps", "trispl")
  				} else {class(rslt) <- c("trispl",  "wps")}	
  			}
  			if (!is.null(warp.delta) & !is.null(delta_add)) {
  				if (labels[1] == "additive") {
  					class(rslt) <- c("cgam", "trispl", "wps")
  				} else if (grepl("warp", labels[1], fixed = TRUE)) {
  					class(rslt) <- c( "wps", "trispl", "cgam")
  				} else {class(rslt) <- c("trispl", "cgam", "wps")} 					
  			}
		}
	} 
	rslt$call <- cl
	rslt$labels <- labels
  	return (rslt) 
}

###############
#amat function#
###############
amat.fun <- function(x)
{
	obs <- 1:length(x)
	amat <- NULL
	xu <- unique(x)
	nx <- sort(xu)
	hd <- head(nx, 1)
	tl <- nx
	while (length(tl) > 1) {
		hd <- head(tl, 1)
		tl <- tl[-1]
		paired <- 0
		for (i in 1:length(tl)) {
			a1 <- 1:length(x)*0
			if (hd * tl[i] > 0) {
				if (hd < 0) {
		        		if (tl[i] > hd) {
						a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1 
					} else {
						a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
					}
				}
				if (hd > 0) {
					if (tl[i] < hd) {
						a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1 
					} else {
						a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
					}
				}
			#if (!all(a1 == 0) & paired == 0 ) {amat <- rbind(amat, a1); paired <- 1}
			}
			if (hd * tl[i] == 0) {
				if (hd == 0) {
					#if (tl[i] > 0) {
						a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
					#}
				} else {
					a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1
				}
			}
			if (!all(a1 == 0) & paired == 0 ) {amat <- rbind(amat, a1); paired <- 1}
		}
	}
	dimnames(amat) <- NULL
	amat
}	

###############
#bmat function#
###############
bmat.fun <- function(x)
{
	obs <- 1:length(x)
	bmat <- NULL
	hd <- head(x, 1)
	tl <- x
	j <- 0
	while (length(tl) > 1) {
		hd <- head(tl, 1)
		tl <- tl[-1]
		paired <- 0
		j <- j + 1
		for (i in 1:length(tl)) {
			b1 <- 1:length(x)*0
			if (hd == tl[i]) {b1[j] <- -1; b1[j + i] <- 1}
			if (!all(b1 == 0) & paired == 0 ) {bmat <- rbind(bmat, b1); paired <- 1}
		}
	}
	dimnames(bmat) <- NULL
	bmat
}

##########
#cgam.fit#
##########
cgam.fit <- function(y, xmat, zmat, shapes, numknots, knots, space, nsim, family = gaussian(), cpar = 1.2, wt.iter = FALSE, umbrella.delta = NULL, tree.delta = NULL, weights = NULL, zid = zid, zid1 = zid1, zid2 = zid2, sc_x = FALSE, sc_y = FALSE, idx_s = NULL, idx = NULL, var.est = NULL, data = NULL) {
  	cicfamily <- CicFamily(family)
	llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
	linkfun <- cicfamily$linkfun
	etahat.fun <- cicfamily$etahat.fun
	gr.fun <- cicfamily$gr.fun
	wt.fun <- cicfamily$wt.fun
	zvec.fun <- cicfamily$zvec.fun
	muhat.fun <- cicfamily$muhat.fun
	ysim.fun <- cicfamily$ysim.fun
	deriv.fun <- cicfamily$deriv.fun
	dev.fun <- cicfamily$dev.fun 
	n <- length(y)
	sm <- 1e-7 
	#sm <- 1e-5
	capl <- length(xmat) / n
	if (capl < 1) {capl <- 0}
	if (round(capl, 8) != round(capl, 1)) {stop ("Incompatible dimensions for xmat!")}
#new:
	if (capl > 0 & sc_x) {
		for (i in 1:capl) {xmat[,i] <- (xmat[,i] - min(xmat[,i])) / (max(xmat[,i]) - min(xmat[,i]))}
		#for (i in 1:capl) {xmat[,i] <- (xmat[,i] - mean(xmat[,i])) / sd(xmat[,i])}
		#for (i in 1:capl) {xmat[,i] <- xmat[,i] / sd(xmat[,i])}
	}
#new:
	if (sc_y) {
		sc <- sd(y)
		y <- y / sc	
	}
	capk <- length(zmat) / n
	if (capk < 1) {capk <- 0}
	if (round(capk, 8) != round(capk, 1)) {stop ("Incompatible dimensions for zmat!")}
#new:
	capls <- sum(shapes == 17)
####################################################
#get basis functions for the constrained components#
####################################################	
	delta <- NULL
	varlist <- NULL
	xid1 <- NULL; xid2 <- NULL; xpos2 <- 0
	knotsuse <- list(); numknotsuse <- NULL
	mslst <- list()
#new:
	capm <- 0
	capms <- 0
	if (capl - capls > 0) {
		del1_ans <- makedelta(xmat[, 1], shapes[1], numknots[1], knots[[1]], space = space[1])
		del1 <- del1_ans$amat
		knotsuse[[1]] <- del1_ans$knots 
		mslst[[1]] <- del1_ans$ms
		numknotsuse <- c(numknotsuse, length(del1_ans$knots))
        m1 <- length(del1) / n
#new code: record the number of columns of del1 if shapes0[1] == 17:
		if (shapes[1] == 17) {capms <- capms + m1}
        var1 <- 1:m1*0 + 1
		xpos1 <- xpos2 + 1
		xpos2 <- xpos2 + m1
		xid1 <- c(xid1, xpos1)
		xid2 <- c(xid2, xpos2)
		if (capl == 1) {
        	delta <- del1
         	varlist <- var1
        } else {
	      	for (i in 2:capl) {
#new code:
	        	del2_ans <- makedelta(xmat[,i], shapes[i], numknots[i], knots[[i]], space = space[i])
				del2 <- del2_ans$amat
				knotsuse[[i]] <- del2_ans$knots
				mslst[[i]] <- del2_ans$ms
				numknotsuse <- c(numknotsuse, length(del2_ans$knots))
				m2 <- length(del2) / n
#new code: record the number of columns of del2 if shapes0[i] == 17:
				if (shapes[i] == 17) {capms <- capms + m2}
				xpos1 <- xpos2 + 1
				xpos2 <- xpos2 + m2
				xid1 <- c(xid1, xpos1)
				xid2 <- c(xid2, xpos2)
				delta <- rbind(del1, del2)
				varlist <- 1:(m1 + m2)*0
				varlist[1:m1] <- var1
				varlist[(m1 + 1):(m1 + m2)] <- (1:m2)*0 + i
				var1 <- varlist
				m1 <- m1 + m2
				del1 <- delta
	      	}
	    }
		if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0 & capk > 0) {
			bigmat <- rbind(1:n*0 + 1, t(zmat), t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
			np <- 1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13)  + capms
		} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0 & capk == 0) {
			bigmat <- rbind(1:n*0 + 1, t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
			np <- 1 + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) + capms
		} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0 & capk > 0) {
			bigmat <- rbind(1:n*0 + 1, t(zmat), delta)
			np <- 1 + capk + capms
		} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0 & capk == 0) {
			bigmat <- rbind(1:n*0 + 1, delta)
			np <- 1 + capms
		} else {
			print ("error in capk, shapes!")
		} 
#new:
		capm <- length(delta) / n - capms
	} else {
	  	if (capk + capls > 0) {
#new:
			if (capls  < 1 & capk > 0) {
          		bigmat <- rbind(1:n*0 + 1, t(zmat))
          		np <- 1 + capk
			} else if (capls > 0) {
				delta <- NULL; varlist <- NULL
				del1_ans <- makedelta(xmat[,1], 17, numknots[1], knots[[1]], space = space[1])
				del1 <- del1_ans$amat
				knotsuse[[1]] <- del1_ans$knots 
				mslst[[1]] <- del1_ans$ms
				numknotsuse <- c(numknotsuse, length(del1_ans$knots))
				m1 <- length(del1) / n
				var1 <- 1:m1*0 + 1
				xpos1 <- xpos2 + 1
				xpos2 <- xpos2 + m1
				xid1 <- c(xid1, xpos1)
				xid2 <- c(xid2, xpos2)
				if (capls == 1) {
        			delta <- del1
         			varlist <- var1
          		} else {
					for (i in 2:capls) {
	        			del2_ans <- makedelta(xmat[,i], 17, numknots[i], knots[[i]], space = space[i])
						del2 <- del2_ans$amat
						knotsuse[[i]] <- del2_ans$knots
						mslst[[i]] <- del2_ans$ms
						numknotsuse <- c(numknotsuse, length(del2_ans$knots))
						m2 <- length(del2) / n
						xpos1 <- xpos2 + 1
						xpos2 <- xpos2 + m2
						xid1 <- c(xid1, xpos1)
						xid2 <- c(xid2, xpos2)
						delta <- rbind(del1, del2)
						varlist <- 1:(m1 + m2)*0
						varlist[1:m1] <- var1
						varlist[(m1 + 1):(m1 + m2)] <- (1:m2)*0 + i
						var1 <- varlist
						m1 <- m1 + m2
						del1 <- delta
	      			}
				}
				if (capk < 1){
					bigmat <- rbind(1:n*0 + 1, delta)
					capms <- length(delta) / n
					np <- 1 + capms
				} else {
					bigmat <- rbind(1:n*0 + 1, t(zmat), delta)
					capms <- length(delta) / n
					np <- 1 + capk + capms 
				}			
			}
        } else {bigmat <- matrix(1:n*0 + 1, nrow = 1); capm <- 0; capms <- 0; np <- 1}
	}
	if (!is.null(umbrella.delta)) {
		bigmat <- rbind(bigmat, umbrella.delta)
		capu <- length(umbrella.delta) / n
	} else {capu <- 0}
	if (!is.null(tree.delta)) {
		bigmat <- rbind(bigmat, tree.delta)
		capt <- length(tree.delta) / n
	} else {capt <- 0}
	if (!is.null(umbrella.delta) | !is.null(tree.delta)) 
		delta_ut <- rbind(umbrella.delta, tree.delta)
	if (capl + capk + capu + capt > 0) {
#		if (capl + capu + capt > 0) {
#new:
		if (capl - capls + capu + capt > 0) {
#new:initialize cvec
			cvec <- NULL
#new: initialize face
            face <- NULL
			if (wt.iter) {
				etahat <- etahat.fun(n, y, fml = family$family)
				gr <- gr.fun(y, etahat, weights, fml = family$family)  
				wt <- wt.fun(y, etahat, n, weights, fml = family$family)     
				cvec <- wt * etahat - gr
			} else {wt <- wt.fun(y, etahat, n, weights, fml = family$family)}
			  	zvec <- zvec.fun(cvec, wt, y, fml = family$family)
#        		gmat <- t(bigmat %*% sqrt(diag(wt)))
#new: avoid memory allocation error
				gmat <- t(bigmat)
				for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
			  	dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
                zsend <- gmat[, 1:np, drop = FALSE] 
                #ans <- coneB(zvec, t(dsend), zsend)
                ans <- coneB(zvec, dsend, zsend)
                face <- ans$face
			  	etahat <- t(bigmat) %*% ans$coefs
#new: for monotinic variance estimation
                vh <- NULL
                kts.var <- NULL
                if (!is.null(var.est)) {
                    #x.var <- var.est
                    #test more:
                    if (!is.null(data)) {
                        nms <- names(data)
                        nm <- attributes(var.est)$nm
                        if (nm %in% nms) {
                            x.var <- data[[nm]]
                            attributes(x.var) <- attributes(var.est)
                        }
                    } else {x.var <- var.est}
                    db.exp <- attributes(x.var)$db.exp
                    kts.var <- attributes(x.var)$var.knots
                    shape.var <- attributes(x.var)$shape
                    iter <- 0
                    diff.var <- 100
                    oldeta <- etahat
                    evec0 <- y - etahat
                    #bigwt <- 10*max(evec0)
                    bigwt <- 1000
                    while(diff.var > 1e-8 & iter < 10) {
                        iter <- iter + 1
                        evec <- y - etahat
                        fit.var <- varest(evec, x.var, shape=shape.var, var.knots=kts.var, db.exp=db.exp)
                        v1 <- fit.var$vhat
                        wt0 <- 1/v1
                        wt0[wt0 > bigwt] <- bigwt
                        wt <- wt.fun(y, etahat, n, weights = wt0, fml = family$family)
                        zvec <- zvec.fun(cvec, wt, y, fml = family$family)
                        gmat <- t(bigmat)
                        for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
                        dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
                        zsend <- gmat[, 1:np, drop = FALSE]
                        ans <- coneB(zvec, dsend, zsend)
                        etahat <- t(bigmat) %*% ans$coefs
                        diff.var <- mean((etahat - oldeta)^2)
                        oldeta <- etahat
                    }
                    kts.var <- fit.var$var.knots
                    vh <- v1
                }
			  	if (wt.iter) {
					muhat <- muhat.fun(etahat, fml = family$family)
			  		diff <- 1
					if (family$family == "binomial") {
						mdiff <- abs(max(muhat) - 1) > sm	
					} else {mdiff <- TRUE}
					nrep <- 0
##########
#iterate!#	
##########
			  		while (diff > sm & mdiff & nrep < n^2){
						oldmu <- muhat	
						nrep <- nrep + 1
						gr <- gr.fun(y, etahat, weights, fml = family$family)	
						wt <- wt.fun(y, etahat, n, weights, fml = family$family) 
						cvec <- wt * etahat - gr
						#zvec <- cvec / sqrt(wt)
						zvec <- zvec.fun(cvec, wt, y, fml = family$family)						
#						gmat <- t(bigmat %*% sqrt(diag(wt)))
						gmat <- t(bigmat)
						for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
						dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
        	  			zsend <- gmat[, 1:np, drop = FALSE]
                        #ans <- coneB(zvec, t(dsend), zsend)
                        ans <- coneB(zvec, dsend, zsend, face = face)
						etahat <- t(bigmat) %*% ans$coefs
						muhat <- muhat.fun(etahat, fml = family$family)
						diff <- mean((muhat - oldmu)^2)	
						mdiff <- abs(max(muhat) - 1)
						if (family$family == "binomial") {
							mdiff <- abs(max(muhat) - 1) > sm	
						} else {mdiff <- TRUE}
					}
			 	}
				yhat <- ans$yhat
				coefskeep <- ans$coefs 
########################
#if capk >= 0, we have:#
########################
				zcoefs <- coefskeep[1:(capk + 1)]
######################
#we will always have:#
######################
				vcoefs <- coefskeep[1:np]
#######################
#if capm > 0, we have:#
#######################
				xcoefs <- NULL
				#if (capm > 0) {
				#	xcoefs <- coefskeep[(np + 1):(np + capm)]
				#}
#new:
				if (capl > 0) {
					xcoefs <- coefskeep[(np - capms + 1):(np + capm)]
				}
#######################
#if capu > 0, we have:#
#######################
				ucoefs <- NULL
				if (capu > 0) {
					ucoefs <- coefskeep[(np + 1 + capm):(np + capm + capu)]
				}
#######################
#if capt > 0, we have:#
#######################
				tcoefs <- NULL
				if (capt > 0) {
					tcoefs <- coefskeep[(np + 1 + capm + capu):(np + capm + capu + capt)]
				}
#########################################################		  
#if we have at least one constrained predictor, we have:#
#########################################################
				thvecs <- NULL
				if (capl > 0) {	
#new code:
					dcoefs <- coefskeep[(np - capms + 1):(np + capm)]
					#dcoefs <- coefskeep[(np + 1):(np + capm)]	
#####################################################
#thvecs is f(x), where x has one of the eight shapes#
#####################################################
					thvecs <- matrix(nrow = capl, ncol = n)
	    			ncon <- 1
	    			for (i in 1:capl) {
	    	  			thvecs[i,] <- t(delta[varlist == i,]) %*% dcoefs[varlist == i]
#new:
	    	  			#if (shapes[i] > 2 & shapes[i] < 5) { 
						if (shapes[i] > 2 & shapes[i] < 5 | shapes[i] > 10 & shapes[i] < 13) { 
            		    	ncon <- ncon + 1
            		    	#thvecs[i,] <- thvecs[i,] + zcoefs[ncon] * xmat[,i]
							thvecs[i,] <- thvecs[i,] + vcoefs[capk + ncon] * xmat[,i]
            	  		}
	    			}
				}
#new:order thvecs back
#print (idx_s)
#print (idx)
				#if (!is.null(idx_s)) { 
				if (length(idx_s) > 0) {
					thvecs0 <- thvecs
					thvecs0[idx_s,] <- thvecs[1:length(idx_s), ]
					#if (!is.null(idx)) {
					if (length(idx) > 0) {
						thvecs0[idx,] <- thvecs[(1+length(idx_s)):capl, ]
					}
					thvecs <- thvecs0
				}
				thvecs_ut <- NULL
				if (capu + capt > 0) {
					thvecs_ut <- t(delta_ut) %*% coefskeep[(np + 1 + capm):(np + capm + capu + capt)]
				}
				if (!is.null(thvecs_ut)) {
					thvecs <- rbind(thvecs, t(thvecs_ut))
				}
#new: problem when not gaussian
				if (sc_y) {
					y <- y*sc
					etahat <- etahat*sc
					for (i in 1:nrow(thvecs)) {
						thvecs[i,] <- thvecs[i,] * sc
					}
				}
	  			etakeep <- etahat
				muhatkeep <- muhat.fun(etakeep, fml = family$family)
				wtkeep <- wt
				df_obs <- sum(abs(coefskeep) > 0)
				#llh <- llh.fun(y, muhat, etahat, n, weights, fml = family$family)
				if (family$family == "Gamma") {
					if ((n - cpar * df_obs) <= 0) {
						phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / df_obs
					} else {
						phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / (n - cpar * df_obs)	
					}
				} else {phihat <- NULL}
				#print (phihat)
				llh <- llh.fun(y, muhatkeep, etakeep, phihat, n, weights, fml = family$family)
				if (family$family == "poisson") {		
					mu0 <- mean(y)
					eta0 <- log(mu0)
				} else {mu0 <- NULL}
                #new: get p-values for smooth components
                    pvs <- NULL
                    s.edf <- NULL
                    bstats <- NULL
                    if (is.numeric(shapes)) {
                        #changed to include shape==17
                        if (all(shapes >= 9 & shapes <= 17)) {
                            for (i in 1:capl) {
                                ansi <- cgam.pv(y=y, xmat=xmat, zmat=zmat, shapes=shapes, delta=delta, np=np, capms=capms, numknotsuse=numknotsuse, varlist=varlist, family=family, weights=weights, test_id=i, nsims=100)
                                pvi <- ansi$pv
                                edfi <- 1.5*sum(xcoefs[varlist == i] > 1e-8)
                                bstati <- ansi$bstat
                                pvs <- c(pvs, pvi)
                                s.edf <- c(s.edf, edfi)
                                bstats <- c(bstats, bstati)
                            }
                        }
                    }
                    #new: get p-values for categorical predictors, not for each level
                    pvsz <- NULL
                    z.edf <- NULL
                    fstats <- NULL
                    if (is.null(weights)) {
                        weights <- 1:n*0 + 1
                    }
                    prior.w <- weights
                    if (capk > 0) {
                        sse1 <- sum(prior.w * (y - muhatkeep)^2)
                        ansi <- cgam.pvz(y=y, bigmat=bigmat, df_obs=df_obs, sse1=sse1, np=np, zid=zid, zid1=zid1, zid2=zid2, muhat = muhatkeep, etahat = etakeep, coefskeep = coefskeep, wt.iter=wt.iter, family=family, weights=weights)
                        pvsz <- ansi$pvs
                        z.edf <- ansi$edfs
                        fstats <- ansi$fstats
                    }
                } else if (capk + capls > 0 & capl - capls + capt + capu == 0) {
                    if (is.null(weights)) {
                        weights <- 1:n*0 + 1
                    }
                    prior.w <- weights
                    vmat <- t(bigmat[1:np, , drop = FALSE])
                    if (wt.iter) {
                        nrep <- 0
                        muhat <- mean(y) + 1:n*0
                        etahat <- linkfun(muhat)
                        diff <- 1
                        if (family$family == "binomial") {
                            mdiff <- abs(max(muhat) - 1) > sm
                        } else {mdiff <- TRUE}
                        while (diff > sm & mdiff & nrep < n^2) {
                            nrep <- nrep + 1
                            oldmu <- muhat
                            zhat <- etahat + (y - muhat) * deriv.fun(muhat, fml = family$family)
                            #w <- diag(as.vector(prior.w / deriv.fun(muhat)))
                            #w <- diag(as.vector(prior.w * (deriv.fun(muhat, fml = family$family))^(-1)))
                            #b <- solve(t(vmat) %*% w %*% vmat) %*% t(vmat) %*% w %*% zhat
                            w <- as.vector(prior.w * (deriv.fun(muhat, fml = family$family))^(-1))
                            tvmat <- t(vmat)
                            for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
                            #print (tvmat)
                            b <- solve(tvmat %*% vmat) %*% tvmat %*% zhat
                            etahat <- vmat %*% b
                            muhat <- muhat.fun(etahat, fml = family$family)
                            diff <- mean((muhat - oldmu)^2)
                            mdiff <- abs(max(muhat) - 1)
                            if (family$family == "binomial") {
                                mdiff <- abs(max(muhat) - 1) > sm
                            } else {mdiff <- TRUE}
                        }
                        zcoefs <- b[1:(capk + 1)]
                        #se.beta <-  sqrt(diag(solve(t(vmat) %*% w %*% vmat)))[1:(capk + 1)]
                        se.beta <- sqrt(diag(solve(tvmat %*% vmat)))[1:(capk + 1)]
                        zstat <- zcoefs / se.beta
                        pvals.beta <-  1 - pchisq(zstat^2, df = 1)
#new: gamma only
                        if (family$family == "Gamma") {
                            phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / (n - np)
                        } else {phihat <- NULL}
                    } else {
                        #w <- diag(prior.w)
                        #b <- solve(t(vmat) %*% w %*% vmat) %*% t(vmat) %*% w %*% y
                        w <- prior.w
                        tvmat <- t(vmat)
                        for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
                        b <- solve(tvmat %*% vmat) %*% tvmat %*% y
                        etahat <- vmat %*% b
                        muhat <- muhat.fun(etahat, fml = family$family)
                        sdhat2 <- sum(prior.w * (y - muhat)^2) / (n - np)
                        zcoefs <- b[1:(capk + 1)]
                        se.beta <- sqrt(diag(solve(tvmat %*% vmat) * sdhat2))[1:(capk + 1)]
                        tstat <- zcoefs / se.beta
                        pvals.beta <-  (1 - pt(abs(tstat), df = n - np)) * 2
#new: for var estimation
                        vh <- NULL
                        kts.var <- NULL
                        if (!is.null(var.est)) {
                            x.var <- var.est
                            db.exp <- attributes(x.var)$db.exp
                            kts.var <- attributes(x.var)$var.knots
                            shape.var <- attributes(x.var)$shape
                            iter <- 0
                            diff.var <- 100
                            oldeta <- etahat
                            evec0 <- y - oldeta
                            #bigwt <- 10*max(evec0)
                            bigwt <- 1000
                            while(diff.var > 1e-8 & iter < 10) {
                                iter <- iter + 1
                                evec <- y - etahat
                                fit.var <- varest(evec, x.var, shape=shape.var, var.knots=kts.var, db.exp=db.exp)
                                v1 <- fit.var$vhat
                                #w <- 1/v1
                                w0 <- 1/v1
                                w0[w0 > bigwt] <- bigwt
                                w <- w0
                                tvmat <- t(vmat)
                                for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
                                b <- solve(tvmat %*% vmat) %*% tvmat %*% y
                                etahat <- vmat %*% b
                                diff.var <- mean((etahat - oldeta)^2)
                                oldeta <- etahat
                                #print (diff.var)
                            }
                            kts.var <- fit.var$var.knots
                            vh <- v1
                        }
                    }
#add thvecs if capls > 0:
                thvecs <- NULL
                if (capls > 0) {
                    thvecs <- matrix(nrow = capls, ncol = n)
                    dcoefs <- b[(capk + 2):np]
                    for (i in 1:capls) {
                        thvecs[i,] <- t(delta[varlist == i,]) %*% dcoefs[varlist == i]
                    }
                }
#new:
                if (sc_y) {
                    y <- y*sc
                    etahat <- etahat*sc
                    muhat <- muhat.fun(etahat, fml = family$family)
                    for (i in 1:nrow(thvecs)) {
                        thvecs[i,] = thvecs[i,] * sc
                    }
                }
                llh <- llh.fun(y, muhat, etahat, phihat, n, weights, fml = family$family)
                df_obs <- np
                dfmean <- np
                rslt <- new.env()
                rslt$family <- family
                rslt$wt.iter <- wt.iter
                #rslt$wt <- diag(w)
                rslt$wt <- w
                rslt$bigmat <- bigmat
                rslt$etahat <- etahat
                rslt$muhat <- muhat
                rslt$d0 <- np
                rslt$capm <- 0
                rslt$capms <- capms
                rslt$capk <- capk
                rslt$capu <- capu
                rslt$capt <- capt
                rslt$xid1 <- xid1 + np - capms
                rslt$xid2 <- xid2 + np - capms
                rslt$dfmean <- dfmean
                rslt$edf0 <- dfmean
                #rslt$llh <- llh
                #if (nsim > 0) {
                if (family$family == "binomial") {
                    nobs <- sum(prior.w)
                } else {nobs <- n}
                if ((nobs - np - cpar * (dfmean - np)) <= 0) {
                    rslt$cic <- llh + log(1 + 2 * dfmean / (dfmean - np))
                } else {
                    #new:
                    if (n<=200){cpar=1.5}
                    rslt$cic <- llh + log(1 + 2 * dfmean / (nobs - np - cpar * (dfmean - np)))
                }
                #rslt$cic <- llh + log(1 + 2 * dfmean / (n - np - 1.5 * (dfmean - np)))
                #}
                rslt$zcoefs <- zcoefs
                rslt$coefs <- b
                rslt$vcoefs <- b
                rslt$xcoefs <- b[(capk + 2):np]
                rslt$se.beta <- se.beta
                rslt$pvals.beta <- pvals.beta
                rslt$dev <- dev.fun(y, muhat, etahat, weights, fml = family$family)$dev
                rslt$dev.null <- dev.fun(y, muhat, etahat, weights, fml = family$family)$dev.null
                rslt$df <- n - np
                rslt$df.null <- n - 1
                #rslt$resid_df_obs <- n - np - 1.5 * (df_obs - np)
                rslt$resid_df_obs <- n - cpar * df_obs
                rslt$vhat <- etahat
                rslt$vmat <- vmat
                rslt$etacomps <- thvecs
                rslt$knots <- knotsuse
                rslt$numknots <- numknotsuse
                rslt$ms <- mslst
                #new: used for predict when there are only shape == 17 x's
                #rslt$xmat <- xmat
                rslt$xmat2 <- xmat
                rslt$knots2 <- knotsuse
                rslt$numknots2 <- numknotsuse
                rslt$ms2 <- mslst
                rslt$vh <- vh
                rslt$kts.var <- kts.var
                return (rslt)
            }
##########
#get cic#
##########
		  if (capl - capls + capu + capt > 0 & nsim > 0) {
	  		dfs <- 1:nsim
#new: initialize face
            face <- NULL
	  		for (isim in 1:nsim) {
				#set.seed(123)
	    		ysim <- ysim.fun(n, mu0, fml = family$family)
		  		if (wt.iter) {
					etahat <- etahat.fun(n, ysim, fml = family$family)
					gr <- gr.fun(ysim, etahat, weights, fml = family$family)
					wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)
					cvec <- wt * etahat - gr
				} else {wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)}
					zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
#            		gmat <- t(bigmat %*% sqrt(diag(wt)))
					gmat <- t(bigmat)
					for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
           			dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
            		zsend <- gmat[, 1:np, drop = FALSE]
                    #ans <- try(coneB(zvec, t(dsend), zsend))
                    ans <- try(coneB(zvec, dsend, zsend))
                    face <- ans$face
					if (class(ans) == "try-error") next 
					if (wt.iter) {
						etahat <- t(bigmat) %*% ans$coefs
						muhat <- muhat.fun(etahat, fml = family$family)
						diff <- 1
						if (family$family == "binomial") {
							mdiff <- abs(max(muhat) - 1) > sm	
						} else {mdiff <- TRUE}
##########
#iterate!#
##########
						nrep <- 0 
						while (diff > sm & nrep < n^2 & mdiff > sm) {
							nrep <- nrep + 1
							oldmu <- muhat	
							gr <- gr.fun(ysim, etahat, weights, fml = family$family)
							wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)
							cvec <- wt * etahat - gr
							#zvec <- cvec / sqrt(wt)
							zvec <- zvec.fun(cvec, wt, y, fml = family$family)
#							gmat <- t(bigmat %*% sqrt(diag(wt)))
							gmat <- t(bigmat)
							for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
							dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
							zsend <- gmat[, 1:np, drop = FALSE]
                            #ans <- try(coneB(zvec, t(dsend), zsend))
                            ans <- try(coneB(zvec, dsend, zsend, face=face))
							if (class(ans) == "try-error") next 
							etahat <- t(bigmat) %*% ans$coefs
							muhat <- muhat.fun(etahat, fml = family$family)
							diff <- mean((muhat - oldmu)^2)	
							if (family$family == "binomial") {
								mdiff <- abs(max(muhat) - 1) > sm	
							} else {mdiff <- TRUE}
						}
					}
	    			dfs[isim] <- sum(abs(ans$coefs) > 0)
	  			}
	  			dfmean <- mean(dfs)
#print (dfmean)
		   	} else if (capl - capls + capu + capt > 0 & nsim == 0) {
				dfmean <- NULL
		   	} 
###################################################
#if the user does not give any predictor, we have:#
###################################################
	} else {
		rslt <- new.env()
		rslt$family <- family 
		rslt$wt.iter <- wt.iter
		rslt$muhat <- 1:n*0 + mean(y)
		rslt$etahat <- linkfun(muhat)
		rslt$dfmean <- 1
		print ("No predictor is provided")
		return (rslt)	
	}
	if (capl - capls + capu + capt > 0) {
#new:
#exclude the case we only have z or unrestricted smooth
		xid1 <- xid1 + np - capms
		xid2 <- xid2 + np - capms 
		#xid1 <- xid1 + np
		#xid2 <- xid2 + np
		rslt <- new.env()		
		rslt$phihat <- phihat
		rslt$family <- family 
		rslt$wt.iter <- wt.iter
		rslt$wt <- wtkeep
		rslt$bigmat <- bigmat  
		rslt$etahat <- etakeep
		rslt$muhat <- muhatkeep 
		rslt$d0 <- np
		rslt$capm <- capm
		rslt$capms <- capms 
		rslt$capk <- capk
		rslt$capu <- capu
		rslt$capt <- capt
		rslt$xid1 <- xid1
		rslt$xid2 <- xid2 
		rslt$coefs <- coefskeep
		rslt$vcoefs <- vcoefs
		rslt$xcoefs <- xcoefs
		rslt$zcoefs <- zcoefs
		rslt$ucoefs <- ucoefs
		rslt$tcoefs <- tcoefs 
		rslt$dfmean <- dfmean
		#rslt$llh <- llh
		#if (nsim > 0) {
		if (!is.null(dfmean)) {
#check!
#new: use dfmean - np if n - np - 1.5 * (dfmean - np) < 0
			if ((n - np - cpar * (dfmean - np)) <= 0) {
				rslt$cic <- llh + log(1 + 2 * dfmean / (dfmean - np))
			} else {
                #new:
                if (n<=200){cpar=1.5}
        		rslt$cic <- llh + log(1 + 2 * dfmean / (n - np - cpar * (dfmean - np)))
			}
		} else {rslt$cic <- NULL}
		rslt$edf0 <- dfmean 
		rslt$etacomps <- thvecs
		vmat <- t(bigmat[1:np, , drop = FALSE])
		rslt$vmat <- vmat
		if (is.null(weights)) {
			weights <- 1:n*0 + 1
		}
		prior.w <- weights
		#w <- diag(as.vector(prior.w / deriv.fun(muhatkeep, fml = family$family)))
		w <- as.vector(prior.w / deriv.fun(muhatkeep, fml = family$family))
###############################################
#the case capk = 0 and capk >= 0 are combined:#
###############################################
#debugged: vhat -> muhat.fun(vhat)
		vhat <- vmat %*% vcoefs	
		muvhat <- muhat.fun(vhat, fml = family$family)
#debugged: yhat -> muhatkeep
		sse1 <- sum(prior.w * (y - muhatkeep)^2)
		sse0 <- sum(prior.w * (y - muvhat)^2)
#new: 
		#if ((n - np - cpar * df_obs) <= 0) {
		if ((n - cpar * df_obs) <= 0) {
			#sdhat2 <- sse1 / (df_obs - np)
			sdhat2 <- sse1 / df_obs
		} else {
			#sdhat2 <- sse1 / (n - np - 1.5 * (df_obs - np))
			#sdhat2 <- sse1 / (n - np - cpar * df_obs)
			sdhat2 <- sse1 / (n - cpar * df_obs)
		}
#debugged: vmat -> vmat and duse
#new: coefskeep include zcoefs; bigmat include vmat
		pmat <- vmat
		capbm <- length(bigmat) / n
		bigmat_nv <- bigmat[(np + 1):capbm, , drop = FALSE]	
		coefs_nv <- coefskeep[(np + 1):capbm]	
		duse <- coefs_nv > 1e-8
		if (sum(duse) >= 1) {
			pmat <- cbind(vmat, t(bigmat_nv[duse, , drop = FALSE]))
		}
		if (wt.iter) {
#new:
			tpmat <- t(pmat)
			for (i in 1:n) {tpmat[,i] <- tpmat[,i] * w[i]}
			se2 <- solve(tpmat %*% pmat)
		} else {
#			se2 <- solve(t(pmat) %*% diag(prior.w) %*% pmat) * sdhat2
			tpmat <- t(pmat)
			for (i in 1:n) {tpmat[,i] <- tpmat[,i] * prior.w[i]}
			se2 <- solve(tpmat %*% pmat) * sdhat2
		}		 		
		se.beta <- 1:(capk + 1)*0
		tstat <- 1:(capk + 1)*0
		pvals.beta <- 1:(capk + 1)*0
		rslt$zcoefs <- zcoefs 	
		for (i in 1:(capk + 1)) {
			se.beta[i] <- sqrt(se2[i,i])
			tstat[i] <- zcoefs[i] / se.beta[i]
#new code: n - np - 1.5 * (df_obs - np) must be positive
			#if ((n - np - cpar * df_obs) <= 0) {
			if ((n - cpar * df_obs) <= 0) {
				pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]),  df_obs)) 
				warning ('Effective degrees of freedom is close to the number of observations! Inference about parametric covariates is not reliable!')
			} else {
				#pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]),  n - np - cpar * df_obs)) 
				pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]),  n - cpar * df_obs)) 
			}
		}
		rslt$se.beta <- se.beta
		rslt$pvals.beta <- pvals.beta
		rslt$sse1 <- sse1
		rslt$sse0 <- sse0
		rslt$dev <- dev.fun(y, muhatkeep, etakeep, weights, fml = family$family)$dev
		rslt$dev.null <- dev.fun(y, muhatkeep, etakeep, weights, fml = family$family)$dev.null
		rslt$df <- n - np 
		rslt$df.null <- n - 1
		#rslt$resid_df_obs <- n - np - cpar * df_obs
		rslt$resid_df_obs <- n - cpar * df_obs
		rslt$df_obs <- df_obs
		rslt$vhat <- vhat
		if (length(knotsuse) == 0) {
			knotsuse <- NULL
		}
#new: order back
		knotsuse2 <- knotsuse
		numknotsuse2 <- numknotsuse
		mslst2 <- mslst
		xmat2 <- xmat
		#if (!is.null(idx_s)) {
		if (length(idx_s) > 0) {	
			knotsuse0 <- knotsuse
			numknotsuse0 <- numknotsuse
			mslst0 <- mslst
			knotsuse0[idx_s] <- knotsuse[1:length(idx_s)]
			numknotsuse0[idx_s] <- numknotsuse[1:length(idx_s)]
			mslst0[idx_s] <- mslst[1:length(idx_s)]
			#if (!is.null(idx)) { 
			if (length(idx) > 0) {
				knotsuse0[idx] <- knotsuse[(1+length(idx_s)):capl]
				numknotsuse0[idx] <- numknotsuse[(1+length(idx_s)):capl]
				mslst0[idx] <- mslst[(1+length(idx_s)):capl]
			}
			knotsuse <- knotsuse0
			numknotsuse <- numknotsuse0
			mslst <- mslst0
		}
# the following has shape = 17 at the beginning: knots2 etc 
		rslt$knots2 <- knotsuse2
		rslt$numknots2 <- numknotsuse2
		rslt$ms2 <- mslst2 
		rslt$xmat2 <- xmat2
		rslt$knots <- knotsuse
		rslt$numknots <- numknotsuse 
		rslt$ms <- mslst
		rslt$capms <- capms
		rslt$cpar <- cpar
        rslt$pvs <- pvs
        rslt$s.edf <- s.edf
        rslt$bstats <- bstats
        #new:
        rslt$pvsz <- pvsz
        rslt$z.edf <- z.edf
        rslt$fstats <- fstats
        rslt$vh <- vh
        rslt$kts.var <- kts.var
		#rslt$sdhat2 <- sdhat2
		return (rslt)
	}
}

###########
#CicFamily#
###########
CicFamily <- function(object,...)UseMethod("CicFamily")
CicFamily <- function(object) {
#temp:
#if (object$family == "ordered") {
#	object$family = "gaussian"
#}
#new:
  if (object$family == "gaussian") {
  	linkfun <- function (mu) mu		
  }
  if (object$family == "poisson") {
  	linkfun <- function (mu) log(mu)			
  }
  if (object$family == "binomial") {
  	linkfun <- binomial()$linkfun
  }
  if (object$family == "Gamma") {
  	#linkfun <- function (mu) log(mu)		
     linkfun <- Gamma(link="log")$linkfun
  }	
  llh.fun <- function(y, muhat = NULL, etahat = NULL, phihat = NULL, n = NULL, weights = NULL, fml = object$family){
    sm <- 1e-7
    #sm <- 1e-5
    if (is.null(weights)) {
	weights <- 1:n*0 + 1
    }
    w <- weights
#new: avoid Inf
    if (fml == "poisson") {
      llh <- 2 * sum(w[w!=0] * (muhat[w!=0] - y[w!=0] * etahat[w!=0])) / n
    }
    if (fml == "binomial") {
      llh <- 0
      if (all(0 <= y) & all(y <= 1)) {
        for (i in 1:n) {
          if (muhat[i] > 0 & muhat[i] < 1) {
            llh <- llh + w[i] * (y[i] * log(muhat[i]) + (1 - y[i]) * log(1 - muhat[i])) 
          }
        }
        llh <- (-2/n) * llh
      } else {
          stop ("y values must be 0 <= y <= 1!")
      }
    }
    if (fml == "gaussian") {
      if (all(w == 1)) {
        llh <- log(sum((y - etahat)^2))
      } else {
          llh <- log(sum(w[w!=0] * (y[w!=0] - etahat[w!=0])^2)) - sum(log(w[w!=0])) / n
      }
    }
    if (fml == "Gamma") {
    	vuhat <- 1 / phihat
    	#print (vuhat)
    	#llh <- 2 * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) / n
    	llh <- 2 / n * (vuhat * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) + n * (log(gamma(vuhat)) - vuhat * log(vuhat)) - (vuhat-1) * sum(log(y)))
    	#print (llh)
    }
    llh 
  }
  
  etahat.fun <- function(n, y, fml = object$family){
    if (fml == "poisson") {
      etahat <- 1:n*0 + log(mean(y)) 
    } 
    if (fml == "binomial") {
      etahat <- 1:n*0 
    }
    if (fml == "Gamma") {
    	etahat <- 1:n*0 + log(mean(y)) 
    }
    etahat
  }

  gr.fun <- function(y, etahat = NULL, weights = NULL, fml = object$family){
    n <- length(y)
    if (is.null(weights)) {
      weights <- 1:n*0 + 1
    }
    w <- weights 
    if (fml == "poisson") {
       gr <- w * (exp(etahat) -  y) 
    }
    if (fml == "binomial") {
       if (all(etahat == 0)) { 
         gr <- w * (1/2 - y)
       } else {
	   gr <- 1:n*0
	   for (i in 1:n) {
	     if (etahat[i] > 100) {
		 gr[i] <- w[i] * (1 - y[i]) 
	     } else {gr[i] <- w[i] * (exp(etahat[i]) / (1 + exp(etahat[i])) - y[i])}
           }         
         }
    }
    if (fml == "Gamma") {
    	gr <- w * (1 - y * exp(-etahat))
    }
    gr
  }

  wt.fun <- function(y, etahat = NULL, n = NULL, weights = NULL, fml = object$family){
    if (is.null(weights)) {
	weights <- 1:n*0 + 1
    }
    w <- weights 
    if (fml == "poisson") {
      wt <-  w * exp(etahat)
    }
    if (fml == "binomial") {
      if (all(etahat == 0)){
        #wt <- 1:n*0 + 1/4
	wt <- w * (1:n*0 + 1/4)
      } else {
	  wt <- 1:n*0
          for (i in 1:n) {
            if (etahat[i] > 100) {
              wt[i] <- 0
            } else {
                wt[i] <- w[i] * exp(etahat[i]) / ((1 + exp(etahat[i]))^2)
              }
          }
        }
    }
    if (fml == "gaussian") {
      wt <- w # (1:n*0 + 1) / w 
    }
    if (fml == "Gamma") {
      wt <-  w * y * exp(-etahat)	
    }
    wt <- as.vector(wt)
    wt 
  }

  zvec.fun <- function(cvec = NULL, wt = NULL, y, sm = 1e-7, fml = object$family) {
    n <- length(y)
    if (fml == "gaussian") {
      #zvec <- y
      zvec <- wt^(1/2) * y
    }
    if (fml == "poisson") {	
      #zvec <- cvec / wt
      zvec <- cvec / sqrt(wt) 
    }
    if (fml == "binomial") {
     zvec = 1:n*0
     zvec[wt == 0] <- 1 / sm
     zvec[wt > 0] <- cvec[wt > 0] / sqrt(wt[wt > 0])
    }
    if (fml == "Gamma") {
      zvec <- cvec / sqrt(wt) 
    }
    zvec 
  }

  muhat.fun <- function(etahat, wt = NULL, fml = object$family){
    n <- length(etahat)
    if (fml == "poisson") {
      muhat <- exp(etahat)
    }
    if (fml == "binomial") {
      muhat <- 1:n*0
      id1 = which(etahat > 100)
      id2 = which(etahat <= 100)
      muhat[id1] = 1
      muhat[id2] = exp(etahat[id2]) / (1 + exp(etahat[id2]))
      #muhat[wt == 0] <- 1
      #for (i in 1:n) {
      #  if (etahat[i] > 100) {
      #      muhat[i] <- 1
      #  } else {
      #    muhat[i] <- exp(etahat[i]) / (1 + exp(etahat[i]))
      #  }
      #}
    }
    if (fml == "gaussian") {
      muhat <- etahat
    }
    if (fml == "Gamma") {
      muhat <- exp(etahat)
    }
   muhat 
  }

  ysim.fun <- function(n, mu0 = NULL, fml = object$family, shp0 = NULL) {
    if (fml == "binomial") {
      ysim <- 1:n*0
      ysim[runif(n) < .5] <- 1
    }
    if (fml == "poisson") {
      if (!is.null(mu0)) {
        ysim <- rpois(n, mu0)
      }
    }
    if (fml == "gaussian") {
      ysim <- rnorm(n)
    }
    if (fml == "Gamma") {
      ysim <- rgamma(n, shape=1)
    }
    ysim 
  }

  deriv.fun <- function(muhat, fml = object$family) {
    if (fml == "binomial") {
		deriv <- 1 / (muhat * (1 - muhat))
    }
    if (fml == "poisson") {
		deriv <- 1 / muhat
    }
    if (fml == "gaussian") {
		deriv <- 1
    }
    if (fml == "Gamma") {
    	deriv <- 1 / muhat
    }
   deriv
  }

 dev.fun <- function(y, muhat, etahat, weights, fml = object$family){
  n <- length(y)
  sm <- 1e-7
  #sm <- 1e-5
  if (is.null(weights)) {
	weights <- 1:n*0 + 1
  }
  w <- weights
  vmat <- matrix(1:n*0 + 1, ncol = 1)
  if (fml == "poisson") {
    #dev <- 2 * sum(w * (y * log(y / muhat) - y + muhat))
	dev <- 0
	for (i in 1:n) {
	  if (y[i] == 0) {
            dev <- dev + 2 * w[i] * muhat[i]
          } else {
            dev <- dev + 2 * w[i] * (y[i] * log(y[i] / muhat[i]) - y[i] + muhat[i])
          }
	}
  }
  if (fml == "binomial") {
        dev <- 0
        for (i in 1:n) {
          if (y[i] == 0 & w[i] != 0) {
            dev <- dev + 2 * w[i] * log(w[i] / (w[i] - w[i] * muhat[i]))
          } else if (y[i] == 1 & w[i] != 0) {
              dev <- dev + 2 * w[i] * log(w[i] / (w[i] * muhat[i]))
          } else if (0 < y[i] & y[i] < 1 & w[i] != 0) {
              dev <- dev + 2 * w[i] * y[i] * log(w[i] * y[i] / (w[i] * muhat[i])) + 2 * (w[i] - w[i] * y[i]) * log((w[i] - w[i] * y[i]) / (w[i] - w[i] * muhat[i]))
          } else {
             stop ("y values must be 0 <= y <= 1!")
          }
       }
  }
  if (fml == "gaussian") {
        dev <- sum(w * (y - muhat)^2)
  }
  if (fml == "Gamma") {
    dev <- 0
	for (i in 1:n) {
	  if (y[i] == 0) {
	  		sm <- 1e-4
            dev <- dev + 2 * w[i] * ((sm - muhat[i]) / muhat[i] - log(sm / muhat[i]))
       } else {
            dev <- dev + 2 * w[i] * ((y[i] - muhat[i]) / muhat[i] - log(y[i] / muhat[i]))
       }
	}
  }
###################
#get null deviance#
###################
  if (fml == "binomial" | fml == "poisson" | fml == "Gamma") {
      diff <- 1
      muhat0 <- mean(y) + 1:n*0
      if (fml == "poisson") {
         etahat0 <- log(muhat0)
      } 
      if (fml == "binomial") {
         etahat0 <- log(muhat0 / (1 - muhat0))
      } 		
      if (fml == "Gamma") {
    	etahat0 <- log(muhat0)
      }
      while (diff > sm) {
        oldmu <- muhat0
		zhat <- etahat0 + (y - muhat0) * deriv.fun(muhat0, fml = fml)		
#	wmat <- diag(as.vector(w / deriv.fun(muhat0, fml = fml)))			
#	b <- solve(t(vmat) %*% wmat %*% vmat) %*% t(vmat) %*% wmat %*% zhat
		wm <- as.vector(w / deriv.fun(muhat0, fml = fml))
        tvmat <- t(vmat)
        for (i in 1:n) {tvmat[,i] <- tvmat[,i] * wm[i]}
        b <- solve(tvmat %*% vmat) %*% tvmat %*% zhat
		etahat0 <- vmat %*% b
		muhat0 <- muhat.fun(etahat0, fml = fml)		
		diff <- mean((muhat0 - oldmu)^2)	
     }
     if (fml == "poisson") {
        #dev.null <- 2 * sum(w * (y * log(y / muhat0) - y + muhat0))
		dev.null <- 0
        for (i in 1:n) {
	  		if (y[i] == 0) {
            	dev.null <- dev.null + 2 * w[i] * muhat0[i]
          	} else {
            	dev.null <- dev.null + 2 * w[i] * (y[i] * log(y[i] / muhat0[i]) - y[i] + muhat0[i])
          	}
		}
     }
     if (fml == "binomial") {
        dev.null <- 0
        for (i in 1:n) {
          if (y[i] == 0 & w[i] != 0) {
            dev.null <- dev.null + 2 * w[i] * log(w[i] / (w[i] - w[i] * muhat0[i]))
          } else if (y[i] == 1 & w[i] != 0) {
              dev.null <- dev.null + 2 * w[i] * log(w[i] / (w[i] * muhat0[i]))
          } else if (0 < y[i] & y[i] < 1 & w[i] != 0) {
              dev.null <- dev.null + 2 * w[i] * y[i] * log(w[i] * y[i] / (w[i] * muhat0[i])) + 2 * (w[i] - w[i] * y[i]) * log((w[i] - w[i] * y[i]) / (w[i] - w[i] * muhat0[i]))
          } else {
              stop ("y values must be 0 <= y <= 1!")
	  	  }
        }
     } 
     if (fml == "Gamma") {
      	dev.null <- 0
		for (i in 1:n) {
	  		if (y[i] == 0) {
	  			sm <- 1e-4
            	dev.null <- dev.null + 2 * w[i] * ((sm - muhat0[i]) / muhat0[i] - log(sm / muhat0[i]))
       		} else {
            	dev.null <- dev.null + 2 * w[i] * ((y[i] - muhat0[i]) / muhat0[i] - log(y[i] / muhat0[i]))
       		}
		}
     }
  }
  if (fml == "gaussian") {
#     wmat <- diag(w)
#     b <- solve(t(vmat) %*% wmat %*% vmat) %*% t(vmat) %*% wmat %*% y
     tvmat <- t(vmat)
     for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
     b <- solve(tvmat %*% vmat) %*% tvmat %*% y
     etahat0 <- vmat %*% b
     muhat0 <- muhat.fun(etahat0, fml = fml)	
     dev.null <- sum(w * (y - muhat0)^2)
  }
  rslt <- new.env()
  rslt$dev <- dev
  rslt$dev.null <- dev.null 
  rslt
  }
  ans <- list(llh.fun = llh.fun, etahat.fun = etahat.fun, gr.fun = gr.fun, wt.fun = wt.fun, zvec.fun = zvec.fun, muhat.fun = muhat.fun, ysim.fun = ysim.fun, deriv.fun = deriv.fun, dev.fun = dev.fun, linkfun = linkfun)
  class(ans) <- "CicFamily"
  return (ans)
}

#######################
#eight shape functions#
####################### 
incr <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 1
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

decr <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 2
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x 
} 

conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 3
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 4
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

incr.conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 5
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

decr.conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 6
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

incr.conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 7
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

decr.conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 8
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

#s.incr <- function(x, numknots = 0, knots = 0, space = "E")
#{
#    cl <- match.call()
#    pars <- match.call()[-1]
#    attr(x, "nm") <- deparse(pars$x)
#    attr(x, "shape") <- 9
#    attr(x, "numknots") <- numknots
#    attr(x, "knots") <- knots
#    attr(x, "space") <- space
#    attr(x, "categ") <- "additive"
#    #class(x) <- "additive"
#    x
#}

#s.decr <- function(x, numknots = 0, knots = 0, space = "E")
#{
#    cl <- match.call()
#    pars <- match.call()[-1]
#    attr(x, "nm") <- deparse(pars$x)
#    attr(x, "shape") <- 10
#    attr(x, "numknots") <- numknots
#    attr(x, "knots") <- knots
#    attr(x, "space") <- space
#    attr(x, "categ") <- "additive"
#    #class(x) <- "additive"
#    x
#}


s.incr <- function(x, numknots = 0, knots = 0, var.knots = 0, space = "E", db.exp = FALSE)
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 9
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    attr(x, "db.exp") <- db.exp
    attr(x, "var.knots") <- var.knots
    #class(x) <- "additive"
    x
}

s.decr <- function(x, numknots = 0, knots = 0, var.knots = 0, space = "E", db.exp = FALSE)
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 10
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    attr(x, "db.exp") <- db.exp
    attr(x, "var.knots") <- var.knots
    #class(x) <- "additive"
    x
}


s.conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 11
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots	
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

s.conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 12
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots	
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

s.incr.conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 13
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots	
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

s.incr.conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 14
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots	
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

s.decr.conv <- function(x, numknots = 0, knots = 0, space = "E") 
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- 15
    attr(x, "numknots") <- numknots
    attr(x, "knots") <- knots	
    attr(x, "space") <- space
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x
}

s.decr.conc <- function(x, numknots = 0, knots = 0, space = "E") 
{
   cl <- match.call()
   pars <- match.call()[-1]
   attr(x, "nm") <- deparse(pars$x)
   attr(x, "shape") <- 16
   attr(x, "numknots") <- numknots
   attr(x, "knots") <- knots	
   attr(x, "space") <- space
   attr(x, "categ") <- "additive"
   #class(x) <- "additive"
   x
}

s <- function(x, numknots = 0, knots = 0, space = "E") 
{
   cl <- match.call()
   pars <- match.call()[-1]
   attr(x, "nm") <- deparse(pars$x)
   attr(x, "shape") <- 17
   attr(x, "numknots") <- numknots
   attr(x, "knots") <- knots	
   attr(x, "space") <- space
   attr(x, "categ") <- "additive"
   #class(x) <- "additive"
   x
}

######################################################
#tree function: give the tree shape to x and return x#
######################################################
tree <- function(x, pl = NULL)
{
    cl <- match.call()
    pars <- match.call()[-1]
    attr(x, "nm") <- deparse(pars$x)
    attr(x, "shape") <- "tree"
    #stop (print (x))
    #print (class(x))
    if (is.null(pl)) {
		if (is.numeric(x)) {
			if (0%in%x) {
				pl <- 0
			} else {pl <- min(x)}
		} else {
			xu <- unique(x)
			pl <- xu[1]
		}
    } else {
		if (!(pl%in%x)) {
			stop ("placebo level is not a level of the tree variable!")
		}
    }
    attr(x, "pl") <- pl
    attr(x, "categ") <- "additive"
    #class(x) <- "additive"
    x

}

##################################################
#tree.fun: make delta to a tree ordering variable#
##################################################
#tree.fun <- function(x)
#{
#    if (min(x) != 0) {
#	stop ("A tree ordering variable must have its placebo equal to 0!")
#    }
#    if (!all(round(x, 0) == x)) {
#	stop ("All elements of a tree ordering variable must be integers!")
#    }
#    if (any(x < 0))
#	stop ("All elements of a tree ordering variable must be positive!")
#    nx <- x
#    obs <- 1:length(x)
#    delta <- matrix(0, nrow = length(attributes(factor(x))$levels) - 1, ncol = length(x))
#    pl <- min(nx)
#    for (i in 1:nrow(delta)) {
#       nx <- nx[which(nx != pl)]
#       pl <- min(nx)
#       index <- obs[x == pl]
#       delta[i, index] <- 1
#    }
#  attr(delta, "shape") <- "tree"
#  delta
#}

tree.fun <- function(x, pl = NULL)
{
  #if (pl %in% x) {
  if (is.null(pl)) {
	if (is.numeric(x)) {
		if (0%in%x) {
			pl <- 0
		} else {pl <- min(x)}
	} else {
		xu <- unique(x)
		pl <- xu[1]
	}
  } else {
	if (!(pl%in%x)) {
		stop ("placebo level is not a level of the tree variable!")
	}
  }
  pos <- x %in% pl
  xu <- unique(x)
  nx <- xu[xu != pl] 
  delta <- matrix(0, nrow = (length(xu) - 1), ncol = length(x))
  for (i in 1:nrow(delta)) {
  	xi = nx[i]
	delta[i, !pos & x %in% xi] <- 1
	pos <- pos | x %in% xi
  }
  #} else {stop ("placebo level is not a level of the tree variable!")}
  attr(delta, "shape") <- "tree"
  delta
}

##############################################################
#umbrella function: give the umbrella shape to x and return x#
##############################################################
umbrella <- function(x)
{
  cl <- match.call()
  pars <- match.call()[-1]
  attr(x, "nm") <- deparse(pars$x)
  attr(x, "shape") <- "umbrella"
  attr(x, "categ") <- "additive"
  #class(x) <- "additive"
  x
}

###########################################################
#umbrella.fun: make delta to an umbrella ordering variable#
###########################################################
umbrella.fun <- function(x)
{
	amat <- amat.fun(x)
	bmat <- bmat.fun(x)
	constr <- t(rbind(amat, bmat))
	vmat <- qr.Q(qr(constr), complete = TRUE)[, -(1:(qr(constr)$rank)), drop = FALSE]
	if (!is.null(bmat)) {
		wperp <- t(rbind(t(vmat), bmat))
		wmat <- qr.Q(qr(wperp), complete = TRUE)[, -(1:(qr(wperp)$rank)), drop = FALSE]	
		atil <- amat %*% wmat 
		delta <- t(wmat %*% t(atil) %*% solve(atil %*% t(atil)))
	} else {
		delta <- t(t(amat) %*% solve(amat %*% t(amat)))
	}
	attr(delta, "shape") <- "umbrella"
	delta
}

###################################################
#find delta for a specific predictor x and a shape#
################################################### 

###################################################
#find delta for a specific predictor x and a shape#
################################################### 
makedelta = function(x, sh, numknots = 0, knots = 0, space = "E", suppre = FALSE, interp = FALSE) {
#if (!interp) {
#x = (x - min(x)) / (max(x) - min(x))
#}
	n = length(x)
# find unique x values
#round(x,8) will make 0 edge in amat!
	#xu = sort(unique(round(x, 8)))
#new: center and scale to avoid numerical instabillity
	xu = sort(unique(x))
	n1 = length(xu)
	sm = 1e-7
	ms = NULL
#  increasing or decreasing
	if (sh < 3) {
		amat = matrix(0, nrow = n1 - 1, ncol = n)
		for (i in 1: (n1 - 1)) {
			amat[i, x > xu[i]] = 1
		}
		if (sh == 2) {amat = -amat}
		if (!interp) {
			for (i in 1:(n1 - 1)) {
#new: use ms in predict.cgam
				ms = c(ms, mean(amat[i, ]))
				amat[i, ] = amat[i, ] - mean(amat[i, ])
			}
		}
	} else if (sh == 3 | sh == 4) {
#  convex or concave
		amat = matrix(0, nrow = n1 - 2 ,ncol = n)
		#for (i in 1: (n1 - 2)) {
		#	amat[i, x > xu[i]] = x[x > xu[i]] - xu[i]
		#}
		for (i in 1: (n1 - 2)) {
			amat[i, x > xu[i+1]] = x[x > xu[i+1]] - xu[i+1]
		}
		if (sh == 4) {amat = -amat}
		xm = cbind(1:n*0+1,x)
		xpx = solve(t(xm) %*% xm)
		pm = xm %*% xpx %*% t(xm)
#new: use ms in predict.cgam
		if (!interp) {
			ms = amat %*% t(pm)
			#amat = amat - amat %*% t(pm)
			amat = amat - ms
		}
	} else if (sh > 4 & sh < 9) {
		amat = matrix(0, nrow = n1 - 1, ncol = n)
		if (sh == 5) { ### increasing convex
			for (i in 1:(n1 - 1)) {
				amat[i, x > xu[i]] = (x[x > xu[i]] - xu[i]) / (max(x) - xu[i])
			}
			if (!interp) {
				for (i in 1:(n1 - 1)) {
					ms = c(ms, mean(amat[i, ]))
					amat[i,] = amat[i,] - mean(amat[i,])
				}
			}
		} else if (sh == 6) {  ## decreasing convex
			for (i in 1:(n1 - 1)) {
				amat[i, x < xu[i + 1]] = (x[x < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
			}
			if (!interp) {
				for (i in 1:(n1 - 1)) {
					ms = c(ms, mean(amat[i, ]))			
					amat[i,] = amat[i,] - mean(amat[i, ])
				}
			}
#print (ms)
		} else if (sh == 7) { ## increasing concave
			for (i in 1:(n1 - 1)) {
				amat[i, x < xu[i + 1]] = (x[x < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
			}
			if (!interp) {
				for (i in 1:(n1 - 1)) {
					ms = c(ms, mean(amat[i, ]))
					amat[i,] = -amat[i,] + mean(amat[i,])
				}		
			}
		} else if (sh == 8) {## decreasing concave
			for (i in 1:(n1 - 1)) {
				amat[i, x > xu[i]] = (x[x > xu[i]] - xu[i]) / (max(x) - xu[i])
			}
			if (!interp) {
				for (i in 1:(n1 - 1)) {
					ms = c(ms, mean(amat[i, ]))
					amat[i,] = -amat[i,] + mean(amat[i,])
				}
			}
		}
	} else if (sh > 8 & sh < 18) {
        #new: add two knots
		#if (all(knots == 0) & numknots == 0) {
		if (length(knots) < 2 & numknots == 0) {
			if (sh == 9 | sh == 10) {#1 2
                #k = trunc(n1^(1/5)) + 4
                if (n1 <= 50) {
                    k = 5
                } else if (n1>50 && n1<100) {
                    k = 6
                } else if (n1>= 100 && n1<200) {
                    k = 7
                } else {
                    k = trunc(n1^(1/5)) + 6
                }
			} else {
                #k = trunc(n1^(1/7) + 4)
                if (n1 <= 50) {
                    k = 5
                } else if (n1>50 && n1<100) {
                    k = 6
                } else if (n1>= 100 && n1<200) {
                    k = 7
                } else {
                    k = trunc(n1^(1/7)) + 6
                }
            }
			if (space == "Q") {
				t = quantile(xu, probs = seq(0, 1, length = k), names = FALSE)
			}
			if (space == "E") {
				#t = 0:k / k * (max(x) - min(x)) + min(x)
				t = 0:(k-1) / (k-1) * (max(x) - min(x)) + min(x)
			} 
		#} else if (any(knots != 0) & numknots == 0) {
		} else if (length(knots) >= 2 & numknots == 0) {
			t = knots
		#} else if (all(knots == 0) & numknots != 0) {
		} else if (length(knots) < 2 & numknots != 0) {
			if (space == "Q") {
				t = quantile(xu, probs = seq(0, 1, length = numknots), names = FALSE)
			} 
			if (space == "E") {
				#k = numknots
#new: numknots should be the # of all knots
				k = numknots - 1
				#if (sh == 9 | sh == 10) {#1 2
				#	k = trunc(n1^(1/5)) + 4
				#} else {k = trunc(n1^(1/7) + 4)}
				t = 0:k / k * (max(x) - min(x)) + min(x)
			}
		#} else if (any(knots != 0) & numknots != 0) {
		} else if (length(knots) >= 2 & numknots != 0) {
			#t0 = quantile(xu, probs = seq(0, 1, length = numknots), names = FALSE)
			t = knots
			if (!suppre) {
				print("'knots' is used! 'numknots' is not used!")
			}
			#print ("'knots' is used!")
			#if (numknots != length(knots)) {
			#	if (!suppre) {
			#		print("length(knots) is not equal to 'numknots'! 'knots' is used!")
			#	}
			#} else if (any(t0 != knots)) {
			#	if (!suppre) {
			#		print("equal x-quantiles knots != 'knots'! 'knots' is used! ") 
			#	}
			#}
		}
		if (sh == 9) {#1			
			amat_ans = monincr(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 10) {#2
			amat_ans = mondecr(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 11) {#3
			amat_ans = convex(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 12) {#4
			amat_ans = concave(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 13) {#5
			amat_ans = incconvex(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 14) {#6
			amat_ans = incconcave(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
		} else if (sh == 15) {#7
			#amat_ans = -incconcave(x, t, interp)
			amat_ans = incconcave(x, t, interp)
			amat = -amat_ans$sigma
			if (!interp) {
				ms = -amat_ans$ms
			}
		} else if (sh == 16) {#8
			#amat_ans = -incconvex(x, t, interp)
			amat_ans = incconvex(x, t, interp)
			amat = -amat_ans$sigma
			if (!interp) {
				ms = -amat_ans$ms
			}
		} else if (sh == 17) {#unconstrained
			amat_ans = incconvex(x, t, interp)
			amat = amat_ans$sigma
			ms = amat_ans$ms
			#amat = -incconcave(x, t)
			#amat = rbind(x, t(bcspl(x, m = length(t), knots = t)$bmat)) 
			#amat = rbind(x, convex(x, t))
		}
	}
	#if (sh < 9) {
	#	rslt = list(amat = amat, knots = 0, ms = ms)
	#} else {
	#	rslt = list(amat = amat, knots = t, ms = ms)
	#}
	if (sh < 9) {t = 0}	
	rslt = list(amat = amat, knots = t, ms = ms)
	rslt
}

#####
#make basis for additive components
make_delta_add = function(xmat, shapes, numknots, knots, space) {
	capl = ncol(xmat) 
	if (capl < 1) {capl = 0}
	if (round(capl, 8) != round(capl, 1)) {stop ("Incompatible dimensions for xmat!")}
	capls = sum(shapes == 17)
	delta = NULL
	varlist = NULL
	knotsuse = list(); numknotsuse = NULL
	mslst = list()
	capm = 0
	capms = 0
	del1_ans = makedelta(xmat[, 1], shapes[1], numknots[1], knots[[1]], space = space[1], interp = FALSE)
	del1 = del1_ans$amat
	knotsuse[[1]] = del1_ans$knots 
	mslst[[1]] = del1_ans$ms
	numknotsuse = c(numknotsuse, length(del1_ans$knots))
    m1 = nrow(del1)
#new code: record the number of columns of del1 if shapes0[1] == 17:
	if (shapes[1] == 17) {capms = capms + m1}
    var1 = 1:m1*0 + 1
	if (capl == 1) {
        delta = del1
        varlist = var1
    } else {
	    for (i in 2:capl) {
	        del2_ans = makedelta(xmat[,i], shapes[i], numknots[i], knots[[i]], space = space[i], interp = FALSE)
			del2 = del2_ans$amat
			knotsuse[[i]] = del2_ans$knots
			mslst[[i]] = del2_ans$ms
			numknotsuse = c(numknotsuse, length(del2_ans$knots))
			m2 = nrow(del2)
#new code: record the number of columns of del2 if shapes0[i] == 17:
			if (shapes[i] == 17) {capms = capms + m2}
			delta = rbind(del1, del2)
			varlist = 1:(m1 + m2)*0
			varlist[1:m1] = var1
			varlist[(m1 + 1):(m1 + m2)] = (1:m2)*0 + i
			var1 = varlist
			m1 = m1 + m2
			del1 = delta
	    }
	}
	if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0) {
		bigmat = rbind(t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
		np = sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13)  + capms
	} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0) {
		bigmat = delta
		np = capms
	} 
	rslt = list(bigmat=bigmat, np=np, varlist=varlist)
	return(rslt)
}

# Monotone increasing
monincr = function(xs, t, interp = FALSE) {
	n = length(xs)
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	k = length(t) - 2
	m = k + 2
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	knt = 1:m
	for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	for (j in 1:(k-1)) {
		index = x >= t[1] & x <= t[j]
		sigma[j, index] = 0

		index = x > t[j] & x <= t[j+1]
		sigma[j, index] = (x[index] - t[j])^2 / (t[j+2] - t[j]) / (t[j+1] - t[j])

		index = x > t[j+1] & x <= t[j+2]
		sigma[j, index] = 1 - (x[index] - t[j+2])^2 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
	    
		index = x > t[j+2] #& x <= t[m]
		sigma[j, index] = 1
	}
	index = x >= t[1] & x <= t[k]
	sigma[k, index] = 0
	
	index = x > t[k] & x <= t[k+1]
	sigma[k, index] = (x[index] - t[k])^2 / (t[k+2] - t[k]) / (t[k+1] - t[k])
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k, index] = 1 - (x[index] - t[k+2])^2 / (t[k+2] - t[k+1]) / (t[k+2] - t[k])
	
	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = 1 - (t[2] - x[index])^2 / (t[2] - t[1])^2

	index = x > t[2] 
	sigma[k+1, index] = 1
	
	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = 0
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = (x[index] - t[k+1])^2 / (t[k+2] - t[k+1])^2
	
#new:
	ms = NULL
	if (!interp) {
		ms = apply(sigma, 1, mean)
		for (i in 1:m) {
			sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	}
	rslt = list(sigma = sigma, ms = ms)
	rslt
}	


########################################################
# Monotone decreasing
mondecr = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	n = length(x)
	k = length(t) - 2
	m = k + 2
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	#knt = 1:m
	#for (i in 1:(k + 2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	#t = x[knt]
	for (j in 1:(k - 1)) {
	 	index = x >= t[1] & x <= t[j]
	 	sigma[j, index] = 1

		index = x > t[j] & x <= t[j+1]
	 	sigma[j, index] = 1 - (x[index] - t[j])^2 / (t[j+2] - t[j]) / (t[j+1] - t[j])

	    	index = x > t[j+1] & x <= t[j+2]
	    	sigma[j, index] = (x[index] - t[j+2])^2 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])

	    	index = x > t[j+2] 
	    	sigma[j, index] = 0
	}

	index = x >= t[1] & x <= t[k]
	sigma[k, index] = 1
	
	index = x > t[k] & x <= t[k+1]
	sigma[k, index] = 1 - (x[index] - t[k])^2 / (t[k+2] - t[k]) / (t[k+1] - t[k])

	index = x > t[k+1] & x <= t[k+2]
	sigma[k, index] = (x[index] - t[k+2])^2 / (t[k+2] - t[k+1]) / (t[k+2] - t[k])

	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = (t[2] - x[index])^2 / (t[2] - t[1])^2

	index = x > t[2] 
	sigma[k+1, index] = 0

	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = 1
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = 1 - (x[index] - t[k+1])^2 / (t[k+2] - t[k+1])^2

	ms = NULL
	if (!interp) {
		ms = apply(sigma, 1, mean)
		for (i in 1:m) {
			sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	}
	rslt = list(sigma = sigma, ms = ms)
	rslt
}

########################################################
# Convex
convex = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	n = length(x)
	k = length(t) - 2
	m = k + 2
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	#knt = 1:m
	#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	for (j in 1:(k-1)) {
	 	index = x >= t[1] & x <= t[j]
	 	sigma[j, index] = 0
	 	
	 	index = x > t[j] & x <= t[j+1]
	 	sigma[j, index] = (x[index] - t[j])^3 / (t[j+2] - t[j]) / (t[j+1] - t[j]) / 3
	    
	   	index = x > t[j+1] & x <= t[j+2]
	    	sigma[j, index] = x[index] - t[j+1] - (x[index] - t[j+2])^3 / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) / 3 + (t[j+1] - t[j])^2 / 3 /(t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])

	    	index = x > t[j+2]
	    	sigma[j, index] = (x[index] - t[j+1]) + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
	}
	index = x >= t[1] & x <= t[k]
	sigma[k, index] = 0
	
	index = x > t[k] & x <= t[k+1]
	sigma[k, index] = (x[index] - t[k])^3 / (t[k+2] - t[k]) / (t[k+1] - t[k]) / 3

	index = x > t[k+1] & x <= t[k+2]
	sigma[k, index] = x[index] - t[k+1] - (x[index] - t[k+2])^3 / (t[k+2] - t[k]) / (t[k+2] - t[k+1]) / 3 + (t[k+1] - t[k])^2 / 3 / (t[k+2] -t[k]) - (t[k+2] - t[k+1])^2 / 3 / (t[k+2] - t[k])
	
	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = x[index] - t[1] + (t[2] - x[index])^3 / (t[2] - t[1])^2 / 3 - (t[2] - t[1]) / 3 #-(t[2]-t[1])^3/(t[2]-t[1])^2/3 #
	
	index = x > t[2] 
	sigma[k+1, index] = x[index] - t[1] - (t[2] - t[1]) / 3 #-(t[2]-t[1])^3/(t[2]-t[1])^2/3#
		
	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = 0
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = (x[index] - t[k+1])^3 / (t[k+2] - t[k+1])^2 / 3
	
	ms = NULL
	if (!interp) {
		xm = cbind(1:n*0+1, x)
		pm = xm %*% solve(t(xm) %*% xm) %*% t(xm)
		ms = matrix(0, nrow = nrow(sigma), ncol = ncol(sigma))
		for (i in 1:m) {
			ms[i,] = pm %*% sigma[i,]
			ms[i,] = ms[i, rank(xs)]		
			#rng=max(sigma[i,])
			#sigma[i,]=sigma[i,]/rng
			sigma[i,] = sigma[i,] - pm %*% sigma[i,]
			sigma[i,] = sigma[i, rank(xs)]
		}
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - pm %*% sigma[i,]
			sigma[i,] = sigma[i, rank(xs)]
		}
	}
	rslt = list(sigma = sigma, ms = ms)
	rslt
}


########################################################
# Concave
concave = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	n = length(x)
	k = length(t) - 2
	m = k + 2
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	#knt = 1:m
	#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	#t = x[knt]
	for (j in 1:k) {
	 	index = x >= t[1] & x <= t[j]
	 	sigma[j, index] = x[index] - t[1]
	 	
	 	index = x > t[j] & x <= t[j+1]
	 	sigma[j, index] = t[j] - t[1] + ((t[j+1] - t[j])^3 - (t[j+1] - x[index])^3) / 3 / (t[j+1] - t[j]) / (t[j+2] - t[j]) + (x[index] - t[j]) * (t[j+2] - t[j+1]) / (t[j+2] - t[j])
	    
	        index = x > t[j+1] & x <= t[j+2]
	    	sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + ((t[j+2] - t[j+1])^3 - (t[j+2] - x[index])^3) / 3 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])	
 	   
 	   	index = x > t[j+2]
 	   	sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
	}

	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = -(t[2] - x[index])^3 / 3 / (t[2] - t[1])^2
	
	index = x > t[2] 
	sigma[k+1, index] = 0
	
	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = x[index] - t[1]
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = t[k+1] - t[1] + ((t[k+2] - t[k+1])^2 * (x[index] - t[k+1]) - (x[index] - t[k+1])^3 / 3) / (t[k+2] - t[k+1])^2
	
	ms = NULL
	if (!interp) {
		xm = cbind(1:n*0+1, x)
		pm = xm %*% solve(t(xm) %*% xm) %*% t(xm)
		ms = matrix(0, nrow = nrow(sigma), ncol = ncol(sigma))
		for (i in 1:m) {
			ms[i,] = pm %*% sigma[i,]
			ms[i,] = ms[i, rank(xs)]		
			sigma[i,] = sigma[i,] - pm %*% sigma[i,]
			sigma[i,] = sigma[i, rank(xs)]
		}
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - pm %*% sigma[i,]
			sigma[i,] = sigma[i, rank(xs)]
		}
	}	
	rslt = list(sigma = sigma, ms = ms)
	rslt
}

########################################################
# Increasing and Convex
incconvex = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	n = length(x)
	k = length(t) - 2
	m = k + 3
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	#knt = 1:(k+2)
	#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	for (j in 1:(k-1)) {
	 	index = x >= t[1] & x <= t[j]
	 	sigma[j, index] = 0
	 	
	 	index = x > t[j] & x <= t[j+1]
	 	sigma[j, index] = (x[index] - t[j])^3 / (t[j+2] - t[j]) / (t[j+1] - t[j]) / 3
	    
	    	index = x > t[j+1] & x <= t[j+2]
	    	sigma[j, index] = x[index] - t[j+1] - (x[index] - t[j+2])^3 / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) / 3 + (t[j+1] - t[j])^2 / 3 /(t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
	    
	  	index = x > t[j+2] 
	    	sigma[j, index] = (x[index] - t[j+1]) + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
	}
	index = x >= t[1] & x <= t[k]
	sigma[k, index] = 0
	
	index = x > t[k] & x <= t[k+1]
	sigma[k, index] = (x[index] - t[k])^3 / (t[k+2] - t[k]) / (t[k+1] - t[k]) / 3
	   
	index = x > t[k+1] & x <= t[k+2]
	sigma[k, index] = x[index] - t[k+1] - (x[index] - t[k+2])^3 / (t[k+2] - t[k]) / (t[k+2] - t[k+1]) / 3 + (t[k+1] - t[k])^2 / 3 / (t[k+2] - t[k]) - (t[k+2] - t[k+1])^2 / 3 / (t[k+2] - t[k])
	
	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = x[index] - t[1] + (t[2] - x[index])^3 / (t[2] - t[1])^2 / 3 - (t[2] - t[1]) / 3
	
	index = x > t[2] 
	sigma[k+1, index] = x[index] - t[1] - (t[2] - t[1]) / 3
	
	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = 0
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = (x[index] - t[k+1])^3 / (t[k+2] - t[k+1])^2 / 3
	
	sigma[k+3,] = x

	ms = NULL
	if (!interp) {
		ms = apply(sigma, 1, mean)
		for (i in 1:m) {
			sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	}
	rslt = list(sigma = sigma, ms = ms)
	rslt
}

########################################################
# Increasing and Concave
incconcave = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
	x = sort(xs)
	n = length(x)
	k = length(t) - 2
	m = k + 3
	sigma = matrix(0, nrow = m, ncol = n)
	obs = 1:n
	#knt = 1:(k+2)
	#for(i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
	for (j in 1:k) {
	 	index = x >= t[1] & x <= t[j]
	 	sigma[j, index] = x[index] - t[1]
	 	
	 	index = x > t[j] & x <= t[j+1]
	 	sigma[j, index] = t[j] - t[1] + ((t[j+1] - t[j])^3 - (t[j+1] - x[index])^3) / 3 / (t[j+1] - t[j]) / (t[j+2] - t[j]) + (x[index] - t[j]) * (t[j+2] - t[j+1]) / (t[j+2] - t[j])
	    
	    	index = x > t[j+1] & x <= t[j+2]
	    	sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + ((t[j+2] - t[j+1])^3 - (t[j+2] - x[index])^3) / 3 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
	   
	    	index = x > t[j+2] 
	    	sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + (t[j+2] - t[j+1])^2 / 3 /(t[j+2] - t[j])
	    
	}

	index = x >= t[1] & x <= t[2]
	sigma[k+1, index] = -(t[2] - x[index])^3 / 3 / (t[2] - t[1])^2
	
	index = x > t[2] 
	sigma[k+1, index] = 0
	
	index = x >= t[1] & x <= t[k+1]
	sigma[k+2, index] = x[index] - t[1]
	
	index = x > t[k+1] & x <= t[k+2]
	sigma[k+2, index] = t[k+1] - t[1] + ((t[k+2] - t[k+1])^2 * (x[index] - t[k+1]) - (x[index] - t[k+1])^3 / 3) / (t[k+2] - t[k+1])^2
	
	sigma[k+3, ] = x

	ms = NULL
	if (!interp) {
		ms = apply(sigma, 1, mean)
		for (i in 1:m) {
			sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	} else {
		for (i in 1:m) {
			#sigma[i,] = sigma[i,] - mean(sigma[i,])
			sigma[i,] = sigma[i, rank(xs)]
		} 
	}
	rslt = list(sigma = sigma, ms = ms)
	rslt
}


##############
#summary.cgam#
##############
summary.cgam <- function(object,...) {
    if (!is.null(object$zcoefs)) {
        family <- object$family
        resid_df_obs <- object$resid_df_obs
        wt.iter <- object$wt.iter
        coefs <- object$zcoefs
        se <- object$se.beta
        tval <- coefs / se
        pvalbeta <- object$pvals.beta
        n <- length(coefs)
        sse0 <- object$SSE0
        sse1 <- object$SSE1
        cic <- object$cic
        deviance <- object$deviance
        null_deviance <- object$null_deviance
        df <- object$df
        null_df <- object$null_df
        zid <- object$zid
        #zid1 <- object$zid1 - 1 - length(shapes)
        #zid2 <- object$zid2 - 1 - length(shapes)
        #new: zid1, zid2 just index zmat not bigmat
        zid1 <- object$zid1
        zid2 <- object$zid2
        tms <- object$tms
        is_param <- object$is_param
        is_fac <- object$is_fac
        vals <- object$vals
        pvs <- object$pvs
        s.edf <- object$s.edf
        bstats <- object$bstats
        if (wt.iter) {
            rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "z.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
            rownames(rslt1)[1] <- "(Intercept)"
            if (n > 1) {
                lzid <- length(zid1)
                for (i in 1:lzid) {
                    pos1 <- zid1[i]; pos2 <- zid2[i]
                    for (j in pos1:pos2) {
                        if (!is_param[i]) {
                            rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")
                        } else {
                            rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")
                        }
                    }
                }
            }
            rslt1 <- as.matrix(rslt1)
            #new:
            rslt2 <- NULL
            if (!is.null(pvs)) {
                rslt2 <- data.frame("edf" = round(s.edf, 4), "mixture of Beta" = round(bstats, 4), "p.value" = round(pvs, 4))
                #rownames(rslt2) <- attributes(tms)$term.labels
                #debugged: check more
                if (!is.null(zid)) {
                    rownames(rslt2) <- (attributes(tms)$term.labels)[-(zid-1)]
                } else {
                    rownames(rslt2) <- (attributes(tms)$term.labels)
                }
            }
            
        } else {
            rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "t.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
            rownames(rslt1)[1] <- "(Intercept)"
            if (n > 1) {
                lzid <- length(zid1)
                for (i in 1:lzid) {
                    pos1 <- zid1[i]; pos2 <- zid2[i];
                    for (j in pos1:pos2) {
                        if (!is_param[i]) {
                            rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")
                        } else {
                            rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")
                        }
                    }
                }
            }
            rslt1 <- as.matrix(rslt1)
            rslt2 <- NULL
            if (!is.null(pvs)) {
                rslt2 <- data.frame("edf" = round(s.edf, 4), "mixture of Beta" = round(bstats, 4), "p.value" = round(pvs, 4))
                #debugged: check more
                if (!is.null(zid)) {
                    rownames(rslt2) <- (attributes(tms)$term.labels)[-(zid-1)]
                } else {
                    rownames(rslt2) <- (attributes(tms)$term.labels)
                }
            }
        }
        #if (!is.null(sse0) & !is.null(sse1)) {
        #rslt2 <- cbind(SSE.Linear = sse0, SSE.Full = sse1)
        #new:
        #	rslt2 <- data.frame("SSE.Linear" = sse0, "SSE.Full" = sse1)
        #	rownames(rslt2)[1] <- ""
        #	ans <- list(call = object$call, coefficients = rslt1, residuals = rslt2, zcoefs = coefs, cic = cic, null_deviance = null_deviance, null_df = null_df, deviance = deviance, df = df, resid_df_obs = resid_df_obs, family = family)
        #	class(ans) <- "summary.cgam"
        #	ans
        #} else {
        ans <- list(call = object$call, coefficients = rslt1, coefficients2 = rslt2, zcoefs = coefs, cic = cic, null_deviance = null_deviance, null_df = null_df, deviance = deviance, df = df, resid_df_obs = resid_df_obs, family = family)
        class(ans) <- "summary.cgam"
        ans
        #}
    } else {
        ans <- list(zcoefs = object$zcoefs)
        class(ans) <- "summary.cgam"
        ans
    }
}

#############
#summary.wps#
#############
summary.wps <- function(object,...) {
	if (!is.null(object$zcoefs)) {
		coefs <- object$zcoefs
		se <- object$se.beta
		#tval <- object$tz
		pvalbeta <- object$pvals.beta
		tval <- coefs / se
		n <- length(coefs)
		#sse0 <- object$SSE0
		#sse1 <- object$SSE1
		zid <- object$zid
#new: zid1, zid2 just index zmat not bigmat
		zid1 <- object$zid1
		zid2 <- object$zid2
		#tms <- object$tms
		#zmat <- object$zmat
		#is_mat <- object$is_mat
		is_param <- object$is_param
		is_fac <- object$is_fac
		vals <- object$vals
		tms <- object$tms
#new: 
		cic <- object$cic
		rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "t.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
		rownames(rslt1)[1] <- "(Intercept)"
		if (n > 1) {
			lzid <- length(zid1)
			for (i in 1:lzid) {
				pos1 <- zid1[i]; pos2 <- zid2[i]
				for (j in pos1:pos2) {
					if (!is_param[i]) {
						rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")						
					} else {
						rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")					
					}	
				}
			}
		}
		rslt1 <- as.matrix(rslt1)
		#if (!is.null(sse0) & !is.null(sse1)) {
		#	rslt2 <- data.frame("SSE.Linear" = sse0, "SSE.Full" = sse1)
#new:
		#	rownames(rslt2)[1] <- ""
			#rslt2 <- as.matrix(rslt2)
		#	ans <- list(call = object$call, coefficients = rslt1, residuals = rslt2, zcoefs = coefs) 
		#	class(ans) <- "summary.wps"
		#	ans
		#} else {
			ans <- list(call = object$call, coefficients = rslt1, zcoefs = coefs, cic = cic)
			class(ans) <- "summary.wps"
			ans
		#}
	} else {
		ans <- list(zcoefs = object$zcoefs)
		class(ans) <- "summary.wps"
		ans
	}
}


#####################
#print.summary.cgam #
#####################
print.summary.cgam <- function(x,...) {
    if (!is.null(x$zcoefs)) {
        #if (!is.null(x$se.beta)) {
        cat("Call:\n")
        print(x$call)
        cat("\n")
        cat("Coefficients:")
        cat("\n")
        printCoefmat(x$coefficients, P.values = TRUE, has.Pvalue = TRUE)
        cat("\n")
        if (x$family$family == "binomial") {
            cat("(Dispersion parameter for binomial family taken to be 1)", "\n")
        }
        if (x$family$family == "poisson") {
            cat("(Dispersion parameter for poisson family taken to be 1)", "\n")
        }
        if (x$family$family == "gaussian") {
            cat("(Dispersion parameter for gaussian family taken to be ", round(x$deviance/x$df,4),")","\n", sep="")
        }
        cat("\n")
        cat("Null deviance: ",round(x$null_deviance,4), "", "on", x$null_df, "", "degrees of freedom", "\n")
        cat("Residual deviance: ",round(x$deviance,4), "", "on", x$resid_df_obs, "", "observed degrees of freedom", "\n")
        #if (is.null(x$cic)) {
        #	message("CIC value is not available when there is no shape-restricted predictor")
        #} else {message("CIC: ", round(x$cic,4))}
        if (!is.null(x$coefficients2)) {
            cat("\n")
            cat("Approximate significance of smooth terms: \n")
            printCoefmat(x$coefficients2, P.values = TRUE, has.Pvalue = TRUE)
        }
        if (!is.null(x$cic)) {
            cat("CIC: ", round(x$cic,4))
        }
        #if (!is.null(x$residuals)) {
        #	cat("==============================================================", "\n")
        #	cat("Call:\n")
        #	print(x$call)
        #	cat("\n")
        #	printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
        #}
    } else {
        print ("No linear predictor is defined")	
        #print ("Residual degree of freedom is negative!")
    }
}

####################
#print.summary.wps #
####################
print.summary.wps <- function(x,...) {
	if (!is.null(x$zcoefs)) {
	#if (!is.null(x$se.beta)) {
		cat("Call:\n")
		print(x$call)
		cat("\n")
		cat("Coefficients:")
		cat("\n")
		printCoefmat(x$coefficients, P.values = TRUE, has.Pvalue = TRUE)
		#if (!is.null(x$residuals)) {
			#cat("==============================================================", "\n")
		#	cat("+--------------------------------------+\n")
  		#	cat("|         SSE.Linear vs SSE.Full       |\n")
#printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
  		#	cat("+--------------------------------------+\n")
			#cat("Call:\n")
			#print(x$call)
			#cat("\n")
		#	printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
		#}
		cat("\n")
		if (!is.null(x$cic)) {
			cat("CIC: ", round(x$cic,4), "\n", sep = "")
		}
	} else {
		print ("No linear predictor is defined")	
	}
}


##############
#predict.cgam#
##############
predict.cgam = function(object, newData, interval = c("none", "confidence", "prediction"), type = c("response", "link"), level = 0.95, n.mix = 500,...) {
#print (is.data.frame(newData))
#print (newData)
#new: 
	family = object$family
	cicfamily = CicFamily(family)
	muhat.fun = cicfamily$muhat.fun	
	if (!inherits(object, "cgam")) { 
	    warning("calling predict.cgam(<fake-cgam-object>) ...")
    }
	if (missing(newData) || is.null(newData)) {
	#if (missing(newData)) {
		etahat = object$etahat
		muhat = muhat.fun(etahat, fml = family$family)
		#ans = list(fit = muhat, etahat = etahat, newbigmat = object$bigmat)
		ans = list(fit = muhat)
		return (ans) 
	}
	if (!is.data.frame(newData)) {
		#newData = as.data.frame(newData)
		stop ("newData must be a data frame!")	
	}
    #shapes = object$shapes
#new: used for ci
	prior.w = object$prior.w
    vh = object$vh
	y = object$y
#new: only use shapes for x != umb or tree
	shapes = object$shapesx
	np = object$d0; capm = object$capm; capk = object$capk; capt = object$capt; capu = object$capu
#new:
	xid10 = object$xid1; xid20 = object$xid2; 
	uid1 = object$uid1; uid2 = object$uid2; tid1 = object$tid1; tid2 = object$tid2 
#new:
	xmat0 = object$xmat0; knots0 = object$knots0; numknots0 = object$numknots0; sps0 = object$sps0; ms0 = object$ms0
	zmat = object$zmat; umb = object$umb; tr = object$tr
#new:
	ztb = object$ztb; zid1 = object$zid1; zid2 = object$zid2; iz = 1
	bigmat = object$bigmat; umbrella.delta = object$umbrella.delta; tree.delta = object$tree.delta
	coefs = object$coefs; zcoefs = object$zcoefs; vcoefs = object$vcoefs; xcoefs0 = object$xcoefs; ucoefs = object$ucoefs; tcoefs = object$tcoefs
	tt = object$tms
	Terms = delete.response(tt)
#model.frame will re-organize newData in the original order in formula
	m = model.frame(Terms, newData)
#print (m)
	newdata = m
#print (head(newdata))
#new:
	newx0 = NULL; newxv = NULL; newx = NULL; newx_s = NULL; newu = NULL; newt = NULL; newz = NULL; newv = NULL
	#newz = list(); iz = 1;
	rn = nrow(newdata)
#print (rn)
#new:
	newetahat = 0; newmuhat = 0
#new: newx_sbasis
	newxbasis = NULL; newx_sbasis = NULL; newubasis = NULL; newtbasis = NULL; newbigmat = NULL
#######################
#local helper function#
#######################
	my_line = function(xp = NULL, y, x, end, start) {
		slope = NULL
		intercept = NULL
		yp = NULL
		slope = (y[end] - y[start]) / (x[end] - x[start])
		intercept = y[end] - slope * x[end]
		yp = intercept + slope * xp
		ans = new.env()
		ans$slope = slope
		ans$intercept = intercept
		ans$yp = yp
		ans
	}
#get shape attributes and elem out of newdata
	for (i in 1:ncol(newdata)) {
		if (is.null(attributes(newdata[,i])$shape)) {
			if (is.factor(newdata[,i])) {
				lvli = levels(newdata[,i])
				ztbi = levels(ztb[[iz]])
				newdatai = NULL
				if (!any(lvli %in% ztbi)) {
					stop ("new factor level must be among factor levels in the fit!")
				} else {
					id1 = which(ztbi %in% lvli)
					#if (any(id1 > 1)) {
#delete the base level
					klvls = length(ztbi)
					if (klvls > 1) {						
						newimat = matrix(0, nrow = rn, ncol = klvls-1)
					        for (i1 in 1:rn) {
							if (newdata[i1,i] != ztbi[1]) {
								id_col = which(ztbi %in% newdata[i1,i]) - 1
								newimat[i1,id_col] = 1
							}
						}
						#for (i1 in 1: klvls) {
						#	which(levels(newdatai) == ztbi[i1])
						#	for (i2 in 1:rn) {
						#		if (levels(newdatai[i2, i1]) == lvli[i2]) {
						#			if (lvli[i2] != ztbi[1]) {
						#				newimat[i2, i1] = 1
						#			}
						#		}
						#	}
						#}
						#ztb_use = ztbi[-1]
						#nc = length(id1[id1 > 1])
						#newdatai = matrix(0, nrow = rn, ncol = nc)
						#for (ic in 1:nc) {
						#	id2 = which(newdata[,i] == ztb_use[ic])
						#	newdatai[id2, ic] = 1
						#}
						newdatai = newimat
					}
				}
				#if (length(levels(newdata[,i])) > 2) {
				#	klvls = length(levels(newdata[,i]))
				#	vals = as.numeric(levels(newdata[,i]))
				#	newimat = matrix(0, nrow = rn, ncol = klvls - 1)
				#	for (i1 in 1: (klvls - 1)) {
				#		for (i2 in 1:rn) {
				#			if (newdatai[i2] == vals[i1 + 1]) {
				#				newimat[i2, i1] = 1
				#			}
				#		}
				#	}
				#	newdatai = newimat
				#}
			} else {
				newdatai = newdata[,i]
			}
			newz = cbind(newz, newdatai)
			iz = iz + 1
#print (head(newz))
		}
		if (is.numeric(attributes(newdata[,i])$shape)) {
			newx0 = cbind(newx0, newdata[,i])
			if ((attributes(newdata[,i])$shape > 2 & attributes(newdata[,i])$shape < 5) | (attributes(newdata[,i])$shape > 10 & attributes(newdata[,i])$shape < 13)) {
				newxv = cbind(newxv, newdata[,i])
			}
		} 
		if (is.character(attributes(newdata[,i])$shape)) {
       			if (attributes(newdata[,i])$shape == "tree") {
				newt = cbind(newt, newdata[,i])
			}
       			if (attributes(newdata[,i])$shape == "umbrella") {
				newu = cbind(newu, newdata[,i])	
			}
     		}
	}
#print (head(newt))
#print (head(newu))
#print (head(newx0))
#print (head(newxv))
#print (head(newz))
#print (head(newv))
#new: separate x and x_s, move shape 17 to the beginning
	if (!is.null(shapes)) {
		if (any(shapes == 17)) {
			kshapes <- length(shapes)
        	obs <- 1:kshapes
        	idx_s <- obs[which(shapes == 17)]; idx <- obs[which(shapes != 17)]	
			newx1 <- newx0
			shapes0 <- 1:kshapes*0
			newx1[,1:length(idx_s)] <- newx0[,idx_s]
			shapes0[1:length(idx_s)] <- shapes[idx_s]   
			if (length(idx) > 0) {
				newx1[,(1 + length(idx_s)):kshapes] <- newx0[,idx]
				shapes0[(1 + length(idx_s)):kshapes] <- shapes[idx]
    			}
			newx0 <- newx1; shapes <- shapes0
		}
#new:xmat is ordinal, xmat_s is smooth
xmat = xmat_s = NULL
newx = newx_s = NULL
knots = NULL
ms_x = ms = NULL
sh_x = sh = NULL
		if (all(shapes < 9)) {
			newx = newx0
			xid1 = xid10; xid2 = xid20
			xmat = xmat0
#new:
			sh_x = shapes 
			ms_x = ms0
		} else if (all(shapes > 8)) {
			newx_s = newx0
			xid1_s = xid10; xid2_s = xid20
			xmat_s = xmat0
			numknots = numknots0
			knots = knots0
			sps = sps0
			sh = shapes
			ms = ms0
		} else if (any(shapes > 8) & any(shapes < 9)) {
			newx = newx0[, shapes < 9, drop = FALSE]
			xmat = xmat0[, shapes < 9, drop = FALSE]
#new:
			sh_x = shapes[shapes < 9]
			ms_x = ms0[shapes < 9]
			xid1 = xid10[shapes < 9]; xid2 = xid20[shapes < 9]

			newx_s = newx0[, shapes > 8, drop = FALSE]
			xmat_s = xmat0[, shapes > 8, drop = FALSE]
			sh = shapes[shapes > 8]
			ms = ms0[shapes > 8]
			xid1_s = xid10[shapes > 8]; xid2_s = xid20[shapes > 8]
			numknots = numknots0[shapes > 8]
			knots = knots0[shapes > 8]
			sps = sps0[shapes > 8]
		}
	}
	#if (!is.null(shapes) & any(shapes == 17)) {
	if (!is.null(shapes)) {
		vcoefs = vcoefs[1:(1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))]
	} else {vcoefs = vcoefs[1:(1 + capk)]}
	if (capk > 0) {
		vcoefs_nz = vcoefs[-c(2:(1 + capk))]
	} else {vcoefs_nz = vcoefs}
	newv = cbind(1:rn*0 + 1, newz, newxv)
	newv_nz = cbind(1:rn*0 + 1, newxv)
#print (newv_nz)
#print (vcoefs)
	#etahat_v = as.vector(newv %*% vcoefs)
	etahat_v_nz = as.vector(newv_nz %*% vcoefs_nz)
#print (etahat_v_nz)
#new:   
	etahat_s = 1:rn*0; newx_sbasis = NULL; xs_coefs = NULL; var_xs = NULL
	if (!is.null(newx_s)) {
 		ks = ncol(newx_s)
		del = NULL
		for (i in 1:ks) {
			xi = xmat_s[,i]
			nxi = newx_s[,i]
			if (any(nxi > max(xi)) | any(nxi < min(xi))) {
				stop ("No extrapolation is allowed in cgam prediction!")
			}
#new: scale accordingly
#nxi = (nxi - min(xi)) / (max(xi) - min(xi))
			pos1 = xid1_s[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
			pos2 = xid2_s[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
			xs_coefs = c(xs_coefs, xcoefs0[pos1:pos2])
			msi = ms[[i]]
			deli_ans = makedelta(nxi, sh[i], numknots[i], knots[[i]], space = sps[i], suppre = TRUE, interp = TRUE)
			deli = deli_ans$amat
			if (sh[i] > 10 & sh[i] < 13) {
				#x = xmat_s[,i]
				xs = sort(xi)
				ord = order(xi)
				nx = length(xi)
				obs = 1:nx
				nr = nrow(deli)
				nc = length(nxi)
				ms2 = matrix(0, nrow = nr, ncol = nc)
				for (i1 in 1:nc) {
					for (i2 in 1:nr) {
						ms2[i2, i1] = my_line(xp = nxi[i1], y = msi[i2, ][ord], x = xs, end = nx, start = 1)$yp 
					}
				}
				deli = deli - ms2
			} else {
				deli = deli - msi
			}
#new:
			var_xs = c(var_xs, 1:nrow(deli)*0 + i)
			del = rbind(del, deli)
		}
		newx_sbasis = t(del)
		etahat_s = as.vector(newx_sbasis %*% xs_coefs)
	}
	etahat_x = 1:rn*0; newxbasis = NULL; xcoefs = NULL; var_x = NULL
	if (!is.null(newx)) {
		kx = ncol(xmat)
		del = NULL
		for (i in 1:kx) {
			xi = xmat[,i]
#new:
			nxi = newx[,i]
			if (any(nxi > max(xi)) | any(nxi < min(xi))) {
				stop ("No extrapolation is allowed in cgam prediction!")
			}
#new: scale accordingly
#nxi = (nxi - min(xi)) / (max(xi) - min(xi))
			shi = sh_x[i]
			pos1 = xid1[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
			pos2 = xid2[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
			xcoefs = c(xcoefs, xcoefs0[pos1:pos2])
			deli = pred_del(xi, shi, nxi, ms_x[[i]])
#new:
			var_x = c(var_x, 1:nrow(deli)*0 + i)
			del = rbind(del, deli)
		}
#new: 
		newxbasis = t(del)
		etahat_x = as.vector(newxbasis %*% xcoefs)
	}
#new:
	etahat_u = 1:rn*0; newuedge = NULL
	if (!is.null(newu)) {	
		for (j in 1:ncol(umb)) {
			u = umb[,j]
			nu = length(u)
			us = sort(u)
			ord = order(u)
			if (any(newu[,j] > max(u)) | any(newu[,j] < min(u))) {
				stop ("no extrapolation is allowed in cgam!")
			}
			#u_edges = t(umbrella.fun(u))
			pos1 = uid1[j]; pos2 = uid2[j]
			u_edges = t(bigmat[pos1:pos2, , drop = FALSE])
			nuedge = ncol(u_edges)
			newuedge0 = matrix(0, nrow = rn, ncol = nuedge)
			for (i in 1:nuedge) {
				ue = u_edges[,i]; udist = round(diff(ue), 4); s = sign(udist)
				ueord = ue[ord]; udistord = round(diff(ueord), 4); sord = sign(udistord)
				obs = 1:(nu-1)
				posin = obs[sord != 0] + 1
				pos = unique(c(1, posin, nu))
				uk = us[pos]
				uek = ueord[pos]
				npos = length(pos)
				nd = length(newu[,j])
				for (l in 1:(npos-1)) {
					if (any(newu[,j] == uk[npos])) {
						ids = which(newu[,j] == uk[npos])
						newuedge0[ids,i] = uek[npos]
					}
					if (any(newu[,j] >= uk[l]) & any(newu[,j] < uk[l+1])) {
						ids = which(newu[,j] >= uk[l] & newu[,j] < uk[l+1])
						newuedge0[ids,i] = uek[l]
					}
				}
			}
			newuedge = cbind(newuedge, newuedge0) 
		}
		newubasis = newuedge
		etahat_u = as.vector(newubasis %*% ucoefs)
	}
# we don't have a newtbasis
# estimate tree
	etahat_t = 1:rn*0 
#print (newt)
	if (!is.null(newt)) {
#each column of  tru is a "table" of all levels of a tree var
		tru = unique(tr)#; placebo = min(tru)
		#tr_etahat = 0
		for (i in 1: ncol(newt)) {
			pos1 = tid1[i] - np - capm - capu
			pos2 = tid2[i] - np - capm - capu
			newtu = unique(newt[,i])
			#tcoefi = tcoefs[pos1:pos2]
#check!	add 0 to be the coef for placebo	
			tcoefi = c(0, tcoefs[pos1:pos2]) 
			if (!all(newtu %in% tru[,i])) {
				stop ("new tree ordering factor must be among the old tree ordering factors!")
			}
			#placebo = min(tru[,i])
#new:
			placebo = object$pl[i]
			if (any(newtu != placebo)) {	
				#id_cf = which(tru[,i] %in% newtu) - 1 #no coef for placebo 
				#tr_etahat = tr_etahat + tcoefi[id_cf]
				#id_cf = which(newtu > placebo)
				id_cf = which(newtu != placebo)
				t_use = newtu[id_cf]
				nt = length(t_use)
				for (it in 1:nt) {
					etahat_t[newt[,i] == t_use[it]] = etahat_t[newt[,i] == t_use[it]] + tcoefi[tru[,i] == t_use[it]]
				}
			}
		}
		#etahat_t = tr_etahat
	} 
#print (etahat_v)
#print (etahat_s)
#print (etahat_x)
#print (etahat_u)
#print (etahat_t)
	etahat_z = 1:rn*0; zcf = NULL
#print (newt)
#print (newdata)
	if (!is.null(newz)) {
#delete the coef for 1
		zcoefs = zcoefs[-1]
		etahat_z = newz %*% zcoefs
		#ztbu = unique(ztb)#; placebo = min(ztbu)
		#lz = ncol(newz)
		#for (i in 1:lz) {
		#	pos1 = zid1[i]
		#	pos2 = zid2[i]
		#	zcoefi = zcoefs[pos1:pos2]
		#	zlvi = zcoefs[pos1:pos2]
		#	etahat_z[newz[,i] == 1] = etahat_z[newz[,i] == 1] + zcoefi
		#}
	} 
#print (etahat_z)
	etahat_v = etahat_v_nz + etahat_z	
	newetahat = etahat_v + etahat_s + etahat_x + etahat_u + etahat_t
	newmuhat = as.vector(muhat.fun(newetahat, fml = family$family))
	if (!is.null(newt)) {
		newtbasis = t(tree.delta[,1:nrow(newData),drop = FALSE])
	}
#print (newv)
	#newbigmat = t(cbind(newv, newx_sbasis, newxbasis, newubasis, newtbasis))
	#ans = list(v = newv, xbs = newxbasis, x_sbs = newx_sbasis, ubs = newubasis, tbs = newtbasis, bigmat = newbigmat, vcoefs = vcoefs, xcoefs = xcoefs, xs_coefs = xs_coefs, ucoefs = ucoefs, etahat = newetahat, muhat = newmuhat)
	#ans = list(muhat = newmuhat)
	#muhat = newmuhat