R/ShapeChange.R

Defines functions find_chpt make_dmat make_rs make_arp CicFamily CicFamily bcsplfirderi bcspl bqspl sl_fun jp.pts ip.kts mode.kts find_pen predict.ShapeChange print.summary.ShapeChange plot.ShapeChange residuals.ShapeChange confint.ShapeChange coef.ShapeChange fitted.ShapeChange ip tp jp changept.fit changept

Documented in changept ip jp tp

changept = function(formula, family = gaussian(), data = NULL, k = NULL, knots = NULL, fir = FALSE, q = 3, pnt = FALSE, pen = 0, arp = FALSE, ci = FALSE, nloop = 1e+3, constr = TRUE, param = TRUE, gcv = FALSE, spl = NULL, dmat = NULL, x1 = NULL, xn = 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!")
  	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")
#print ('ev0')
  	mf = eval(mf, parent.frame())
#print ('ev1')
  	ynm = names(mf)[1]
  	mt = attr(mf, "terms")
#ptm=proc.time()
    	y = model.response(mf, "any")
#print (proc.time()-ptm)
  	if (!is.null(names(y))) {
    		y = unname(y)
  	}
	if (family$family == "binomial") {
		#if (class(y) == "factor") {
	  if (inherits(y, "factor")) {
			y = ifelse(y == levels(y)[1], 0, 1)
		}
  	}
	if (family$family == "binomial" | family$family == "poisson") {
		wt.iter = TRUE
	} else {
		wt.iter = FALSE
	}
	zmat = NULL; zid = NULL; zid0 = NULL; zid1 = NULL; zid2 = NULL; znms = NULL; is_fac = NULL; vals = NULL; st = 1; ed = 1
	for (i in 2:ncol(mf)) {
#print (i)
		if (is.character(attributes(mf[, i])$categ)) {
  			x = mf[, i]
  			categ = attributes(x)$categ
  			sh = attributes(x)$sh
			trd1 = attributes(x)$trd1
  			trd2 = attributes(x)$trd2
			up = attributes(x)$up
			xnm = attributes(x)$nm
		}
		if (is.null(attributes(mf[, i])$categ)) {
			if (!is.null(names(mf)[i])) {
				znms = c(znms, names(mf)[i])
				#vals = c(vals, unique(mf[, i]))
			}
			if (!is.matrix(mf[, i])) {
      				zid = c(zid, i)
				is_fac = c(is_fac, 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)
    			  	} else {
					#is_fac = c(is_fac, TRUE)
					zfac = factor(mf[, i])
					ch_char = suppressWarnings(is.na(as.numeric(levels(zfac))))
	        			if (any(ch_char)) {
	    					vals = c(vals, unique(levels(zfac))[-1])
          				} else {
	    					vals = c(vals, as.numeric(levels(zfac))[-1])
	  				}
					zmat = cbind(zmat, model.matrix(~ zfac)[, -1, drop = FALSE])	
					nlvs = length(attributes(zfac)$levels)
	  				ed = st + nlvs - 2 
	  				zid1 = c(zid1, st)
	  				zid2 = c(zid2, ed)
	  				st = st + nlvs - 1
            			}
			} else {
#print ('TRUE')
	 			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)
	  			#vals <- c(vals, 1)
	  			zid = c(zid, i)
#already exclude the base level
	  			nlvs = ncol(zmat0) + 1
#print (nlvs)
	  			ed = st + nlvs - 2
	  			zid1 = c(zid1, st)
	  			zid2 = c(zid2, ed)
	  			st = st + nlvs - 1
			}
		}
	}
	if (is.matrix(zmat) & !is.null(zmat)){
    		nzmat = zmat
    		if (qr(nzmat)$rank != ncol(nzmat)) {
       			stop("zmat should be full column rank!")
    		}
    		mat_cols = ncol(nzmat)
    		mat_rm = NULL
    		for (i in 1:mat_cols) {
       			if (all(round(diff(nzmat[, i]), 8) == 0)) {
              			mat_rm = c(mat_rm, i)
       			}
    		}
    		if (!is.null(mat_rm)) {
           		 nzmat = nzmat[, -mat_rm, drop = FALSE]
    		}
    		zmat = nzmat
  	}
#print ('start0')
  	ans = changept.fit(x, y, zmat = zmat, family = family, categ = categ, m = k, knots = knots, sh = sh, fir = fir, q = q, pnt = pnt, pen = pen, arp = arp, ci = ci, nloop = nloop, trd1 = trd1, trd2 = trd2, up = up, constr = constr, wt.iter = wt.iter, param = param, gcv = gcv, spl = spl, dmat = dmat, x1 = x1, xn = xn)
  	rslt = list(chpt = ans$chpt, knots = ans$knots, fhat = ans$fhat, fhat_x = ans$fhat_x, fhat_eta = ans$fhat_eta, fhat_eta_x = ans$fhat_eta_x, coefs = ans$bhat, zcoefs = ans$zcoefs, cibt = ans$cibt, categ = categ, sh = sh, x = x, y = y, xnm = xnm, znms = znms, ynm = ynm, m_i = ans$m_i, pos = ans$pos, sub = ans$sub, family = family, wt.iter = wt.iter, zmat = zmat, vals = vals, is_fac = is_fac, zid = zid, zid1 = zid1, zid2 = zid2, tms = mt, msbt = ans$msbt, bmat = ans$bmat, phi = ans$phi, sig = ans$sig, aics = ans$aics, lambda = ans$lambda, lam_arp = ans$lams, edfs = ans$edfs, gcvus = ans$gcvus, gcvu = ans$gcvu, lams = ans$lams, phisbt = ans$phisbt, sigsbt = ans$sigsbt, edfcs = ans$edfcs, edfc = ans$edfc, edfu = ans$edfu, gcvcs = ans$gcvcs, ms_gcv = ans$ms_gcv, fhats_gcv = ans$fhats_gcv, fhats_x_gcv = ans$fhats_x_gcv, ssegcvs = ans$ssegcvs, id_gcv = ans$id_gcv, spl = ans$spl, dmat = ans$dmat, x1 = ans$x1, xn = ans$xn)
  	rslt$call = cl
  	class(rslt) = "ShapeChange"
  	return (rslt) 
}
########
changept.fit = function(x, y, zmat = zmat, family = gaussian(), categ = categ, m = NULL, knots = NULL, sh = 1, fir = FALSE, q = 3, pnt = FALSE, pen = 0, arp = FALSE, lambdas = NULL, ci = FALSE, nloop = 1e+3, trd1 = -1, trd2 = -1, up = TRUE, constr = TRUE, wt.iter = wt.iter, param = TRUE, hs0 = NULL, ids0 = NULL, gcv = FALSE, spl = NULL, dmat = NULL, x1 = NULL, xn = NULL) {
	n = length(y)
	if (is.null(hs0) && arp) {  	
     		hs0 = 0:(n - 1)
     		id_1 = id_2 = hs0
     		id_mat = sapply(id_1, function(id_1) id_2 - id_1)
     		ids0 = sapply(hs0, function(h) {which(abs(id_mat) == h)})
	}
	if (categ == "inflect") {
		ans = ip.kts(x, y, zmat = zmat, m = m, knots = knots, q = q, pen = pen, pnt = pnt, arp = arp, lambdas = lambdas, pen_bt = NULL, p_bt = NULL, sh = sh, fir = fir, wt.iter = wt.iter, hs = hs0, ids = ids0, family = family, gcv = gcv, spl = spl, dmat = dmat, x1 = x1, xn = xn) 
	} else if (categ == "mode") {
		ans = mode.kts(x, y, zmat = zmat, m = m, knots = knots, q = q, pen = pen, pnt = pnt, arp = arp, lambdas = lambdas, pen_bt = NULL, p_bt = NULL, sh = sh, wt.iter = wt.iter, hs = hs0, ids = ids0, family = family, gcv = gcv, spl = spl, dmat = dmat, x1 = x1, xn = xn) 
	} else if (categ == "jp") {
		minsse = sum(y^2)
		for (i in 1:(n - 1)) {
			ansi = jp.pts(x, y, zmat = zmat, jpt = (sort(x))[i], i = NULL, m = m, knots = knots, q = q, pen = pen, pnt = pnt, up = up, trd1 = trd1, trd2 = trd2, constr = constr, wt.iter = wt.iter, arp = arp, lambdas = lambdas, pen_bt = NULL, p_bt = NULL, hs = hs0, ids = ids0, family = family, gcv = gcv) 		
			ssei = sum((y - ansi$fhat)^2)
			if (ssei < minsse) { 
				minsse = ssei; ans = ansi  
			}
		}
	} else {
		stop("change-point cannot be estimated!")
	}
	cibt = msbt = trs = trrs = lamsbt = sigsbt = psbt = NULL
	phisbt = list()
	if (ci & nloop > 0) {
		msbt = 1:nloop*0
		trs = trrs = lamsbt = sigsbt = psbt = 1:nloop*0
		phisbt = list()
		yvec = ans$fhat
#new
    		ord = order(x)
    		xs = sort(x)
    		x = xs
    		y = y[ord]
		evec = y - yvec
    		ys = matrix(0, nrow = nloop, ncol = n)
		ans_pen = ans_p = ans_lam = NULL
		for (j in 1:nloop) {
			if (family$family == "gaussian") {
           			if (!arp) {
				      esim = sample(evec, size = n, replace = TRUE)
           			} else {
             				ans_pen = ans$lambda
					ans_phi = ans$phi
					ans_sig = ans$sig
             				ans_p = length(ans$phi)
             				ans_lam = ans$lams
             				if (param) {
                				esim = as.vector(arima.sim(n = n, list(order = c(ans_p, 0, 0), ar = ans_phi), sd = sqrt(ans_sig)))
             				} else {
						#r0 = ans$r0
               					rmat = ans$rmat
						#rmat = rmat * r0
               					yvec = ans$fhat
               					evec = matrix(ans$es, ncol = 1)
               					umat = chol(rmat)
               					e_iid = as.vector(solve(t(umat), evec))
               					esamp = matrix(sample(e_iid, size = n, replace = TRUE), ncol = 1)
               					#esim = sqrt(ans_sig) * crossprod(umat, esamp)
             					esim = crossprod(umat, esamp)
						}
           			}
				ysim = yvec + esim
			} else if (family$family == "binomial") {
				ysim  = rbinom(n, size = 1, prob = yvec)
			} else if (family$family == "poisson") {
				ysim  = rpois(n, lambda = yvec)
			} 
			if (categ == "inflect") {
				ansj = try(ip.kts(x, ysim, zmat = zmat, m = m, knots = knots, q = q, pen = pen, pnt = pnt, arp = arp, lambdas = ans_lam, pen_bt = ans_pen, p_bt = ans_p, sh = sh, fir = fir, wt.iter = wt.iter, hs = hs0, ids = ids0, family = family, gcv = gcv, spl = spl, dmat = dmat, x1 = x1, xn = xn))
        			if (inherits(ansj, "try-error")) next 
			} else if (categ == "mode") {            
				ansj = try(mode.kts(x, ysim, zmat = zmat, m = m, knots = knots, q = q, pen = pen, pnt = pnt, arp = arp, lambdas = ans_lam, pen_bt = ans_pen, p_bt = ans_p, sh = sh, wt.iter = wt.iter, hs = hs0, ids = ids0, family = family, gcv = gcv, spl = spl, dmat = dmat, x1 = x1, xn = xn))
        			if (inherits(ansj, "try-error")) next 
			} else if (categ == "jp") {
				minsse = sum(ysim^2)				
				for (i in 1:(n - 1)) {
					ansi = try(jp.pts(x, ysim, zmat = zmat, jpt = sort(x)[i], i = NULL, m = m, knots = knots, q = q, pen = pen, pnt = pnt, up = up, trd1 = trd1, trd2 = trd2, constr = constr, wt.iter = wt.iter, arp = arp, lambdas = lambdas, pen_bt = ans_pen, p_bt = ans_p, hs = hs0, ids = ids0, family = family, gcv = gcv))
					if (inherits(ansj, "try-error")) next 
#new:
					if (arp) {
						rmat = ansi$rmat 
						ssei = crossprod((ysim - ansi$fhat), chol2inv(chol(rmat))) %*% (ysim - ansi$fhat)
					} else {
						ssei = sum((ysim - ansi$fhat)^2)
					}
					if (ssei < minsse) {
						minsse = ssei; ansj = ansi  
					}
				}
			}
			msbt[j] = ansj$chpt
			if (arp) {
				sigsbt[j] = ansj$sig
				phisbt[[j]] = ansj$phi
				lp = length(ansj$phi)
				if (lp > 1) {
					psbt[j] = max(which(round(ansj$phi, 6) != 0)) 
				} else {
					if (round(ansj$phi, 6) != 0) {
						psbt[j] = 1
					} else {psbt[j] = 0}
				}    						
				trs[j] = ansj$tr
				trrs[j] = ansj$trr
      				lamsbt[j] = ansj$lambda
			}
  	}
	cibt = quantile(sort(msbt), probs = c(.025, .975))
	}
	rslt = list(chpt = ans$chpt, knots = ans$knots, fhat = ans$fhat, fhat_x = ans$fhat_x, fhat_eta = ans$fhat_eta, fhat_eta_x = ans$fhat_eta_x, sse = ans$sse, df = ans$df, pos = ans$pos, dist = ans$dist, bhat = ans$bhat, zcoefs = ans$zcoefs, pvals.beta = ans$pvals.beta, se.beta = ans$se.beta, tval = ans$tval, cibt = cibt, m_i = ans$m_i, sub = ans$sub, msbt = msbt, phisbt = phisbt, psbt = psbt, sigsbt = sigsbt, bmat = ans$bmat, dmat = ans$dmat, phi = ans$phi, sig = ans$sig, aics = ans$aics, aiclst = ans$aiclst, es = ans$es, lambda = ans$lambda, lams = ans$lams, edfs = ans$edfs, gcvus = ans$gcvus, gcvu = ans$gcvu, rmat = ans$rmat, aicmat = ans$aicmat, sigmat = ans$sigmat, trsmat = ans$trsmat, trrsmat = ans$trrsmat, psmat = ans$psmat, fsmat = ans$fsmat, phismat = ans$phismat, psmat = ans$psmat, pensmat = ans$pensmat, hs = hs0, ids = ids0, edfcs = ans$edfcs, edfc = ans$edfc, edfu = ans$edfu, gcvcs = ans$gcvcs, ms_gcv = ans$ms_gcv, fhats_gcv = ans$fhats_gcv, fhats_x_gcv = ans$fhats_x_gcv, ssegcvs = ans$ssegcvs, id_gcv = ans$id_gcv, spl = ans$spl, dmat = ans$dmat, x1 = ans$x1, xn = ans$xn)
	rslt
}

######################
jp = function(x, trd1 = -1, trd2 = -1, up = TRUE) {
  cl = match.call()
  pars = match.call()[-1]
  attr(x, "nm") = deparse(pars$x)
  attr(x, "categ") = "jp"
  attr(x, "trd1") = trd1
  attr(x, "trd2") = trd2
  attr(x, "up") = up
  x
}

tp = function(x, sh = 1) {
  cl = match.call()
  pars = match.call()[-1]
  attr(x, "nm") = deparse(pars$x)
  attr(x, "categ") = "mode"
  attr(x, "sh") = sh
  x
}

ip = function(x, sh = 1) {
  cl = match.call()
  pars = match.call()[-1]
  attr(x, "nm") = deparse(pars$x)
  attr(x, "categ") = "inflect"
  attr(x, "sh") = sh
  x
}

######################
fitted.ShapeChange = function(object, ...) {
  if (!inherits(object, "ShapeChange")) { 
    warning("calling fitted.ShapeChange(<fake-ShapeChange-object>) ...")
  }
  ans = as.vector(object$fhat)
  ans
}

coef.ShapeChange = function(object, ...) {
  if (!inherits(object, "ShapeChange")) { 
    warning("calling coef.ShapeChange(<fake-ShapeChange-object>) ...")
  }
  ans = object$coefs
  ans	
}

confint.ShapeChange = function(object, parm = NULL, level = NULL, ...) {
  if (!inherits(object, "ShapeChange")) { 
    warning("calling confint.ShapeChange(<fake-ShapeChange-object>) ...")
  }
  ans = object$cibt
  ans	
}

residuals.ShapeChange = function(object, ...) {
  if (!inherits(object, "ShapeChange")) { 
    warning("calling resid.ShapeChange(<fake-ShapeChange-object>) ...")
  }
  ans = object$y - object$fhat
  ans	
}

##########################
plot.ShapeChange = function(x, y = NULL, tt = TRUE, xlab = NULL, ylab = NULL, color = "mediumorchid4", type = "p", lty = 1, lwd = 2, cex = .5, has.leg = FALSE, loc = NULL, rugged = FALSE, detailed = FALSE,...) {
	if (!inherits(x, "ShapeChange")) { 
		warning("calling plot.ShapeChange(<fake-ShapeChange-object>) ...")
 	}
 	object = x
 	x = object$x
	zmat = object$zmat
	znm = object$znm
	zid = object$zid
	zid1 = object$zid1 
	zid2 = object$zid2
	zbhat = object$zcoefs
  	tms = object$tms
	vals = object$vals
	znm_2 = dimnames(zmat)[[2]]
	xnm = object$xnm
	ord = order(x)
	x0s = sort(x)
	x = (x0s - min(x0s)) / (max(x0s) - min(x0s))
	y = object$y
	ynm = object$ynm
	y = y[ord]
	sub = object$sub
	fhat = object$fhat
	fhat_x = object$fhat_x
	fhat_eta_x = object$fhat_eta_x
	is_fac = object$is_fac
	bhat = object$coefs
  	chpt = object$chpt	
	pos = object$pos
#constr = object$constr
	ci = object$ci
	m_i = object$m_i
  	categ = object$categ
	sh = object$sh
	knots = object$knots
	knots = (max(x0s) - min(x0s)) * knots + min(x0s)
	k = length(knots)
	lb = length(bhat)
	wt.iter = object$wt.iter
	family = object$family
#print (family$family)
	cicfamily = CicFamily(family)
	muhat.fun = cicfamily$muhat.fun
	if (!is.null(zmat)) {
		zu = sort(unique(zmat))
	}
	palette = c("lightblue", "peachpuff", "limegreen", "grey", "wheat", "yellowgreen", "seagreen1", "palegreen", "azure", "whitesmoke")
  	ylim = c(min(c(min(y), min(fhat))), max(c(max(y), max(fhat))))
#new:
	if (categ == "mode" | categ == "jp") {
		ans0 = bqspl(x, k, pic = TRUE, spl = FALSE)
	} else if (categ == "inflect") {
		ans0 = bcspl(x, k, pic = TRUE, spl = FALSE)
	}
	xpl = ans0$xpl
	xpl0 = 0:1000/1000*(max(x0s) - min(x0s)) + min(x0s)
	bpl = ans0$bpl    
	if (categ == "jp") {
		n = length(x)
		ijp = pos
		d1 = d2 = 1:1001*0
		d1[xpl <= (x[ijp] + x[ijp + 1]) / 2] = (2 * ijp + 1) / 2 - n
		d1[xpl > (x[ijp] + x[ijp + 1]) / 2] = (2 * ijp + 1) / 2
		d2[xpl <= (x[ijp] + x[ijp + 1]) / 2] = 0
		d2[xpl > (x[ijp] + x[ijp + 1]) / 2] = xpl[xpl > (x[ijp] + x[ijp+1]) / 2] - (x[ijp] + x[ijp+1]) / 2
		d2 = d2 - m_i
		if (ijp == 1 | ijp == (n - 1)) {
			d2 = NULL
		}
		di = cbind(d1, d2)	
		bpl = cbind(bpl, di)
	} 
	if (is.null(zmat)) {
		thp = bpl %*% bhat
	} else {
		thp = bpl %*% as.vector(bhat)[1:(ncol(bpl))]	
	}	
	if (tt) {
    		if (categ == "mode") {
     	 		tt0 = "Mode Estimate: "
    	} else if (categ == "jp") {
      		tt0 = "Jump-Point Estimate: "
    	} else if  (categ == "inflect") {
      		tt0 = "Inflection-Point Estimate: "
    	} else {
      		stop ("Wrong Shape!!")
    	}
    		main = paste(tt0, round(chpt, 2))
  	} else { 
    		main = NULL
  	}
	thp0 = thp
	if (wt.iter & categ != "inflect") {
		thp = muhat.fun(thp0, fml = family$family)
	}
	r1 = range(y); r2 = range(thp); r3 = range(fhat)
	rmin = min(c(r1, r2, r3)); rmax = max(c(r1, r2, r3))
	rg = rmax - rmin 
  	if (!is.null(xlab)) {
     		xnm = xlab
  	}
  	if (!is.null(ylab)) {
     		ynm = ylab
  	}
	plot(x0s, y, type = "n", col = 1, cex = cex, ylab = ynm, xlab = xnm, ylim = c(rmin, rmax), main = main)
	if (family$family == "binomial" & rugged) {
		#loc = "right"
		rug(x0s[y == 0])
		rug(x0s[y == 1], side = 3)
	} else {
		points(x0s, y, lwd = 1, cex = cex)
	}	
	if (is.null(zmat)) {
		if (detailed) {
			lines(xpl0, thp, type = 'l', lty = 2, xlab = '', ylab = '', col = color, lwd = 2)
			points(x0s, fhat, col = color, cex = cex)   	
		} else {
  			lines(x0s, fhat, type = 'l', lty = 2, xlab = '', ylab = '', col = color, lwd = lwd)
		}		
	} else {
		lzb = length(zbhat)
		if (detailed) {
			lines(xpl0, thp, type = 'l', lty = 2, xlab = '', ylab = '', col = color, lwd = 2)	
			points(x0s, fhat_x, col = color, cex = cex)
			for (i in 1:lzb) {                                      
				thpi = thp + zbhat[i]
				if (wt.iter & categ != "inflect") {
					thpi = muhat.fun(thpi, fml = family$family)
					fhat_eta_xi = fhat_eta_x + zbhat[i]
					fhat_xi = muhat.fun(fhat_eta_xi, fml = family$family)				
				} else {
					fhat_xi = fhat_x + zbhat[i]
				}
				lines(xpl0, thpi, type = 'l', lty = 2, xlab = '', ylab = '', col = palette[i], lwd = 2)	
				points(x0s, fhat_xi, col = palette[i], cex = cex)			
			}
		} else {
  			lines(x0s, fhat_x, type = 'l', lty = 2, col = color, lwd = 2)
			for (i in 1:lzb) {                                      
				if (wt.iter & categ != "inflect") {
					fhat_eta_xi = fhat_eta_x + zbhat[i]
					fhat_xi = muhat.fun(fhat_eta_xi, fml = family$family)				
				} else {
					fhat_xi = fhat_x + zbhat[i]
				}
				lines(x0s, fhat_xi, type = 'l', lty = 2, col = palette[i], lwd = 2)				
			}
		}
	}
	#text(chpt, ylim[1], round(chpt, 4), col = color, lwd = 2)
	if (!is.null(ci)) {
		lvl = 95
		if (is.null(loc)) {
			if (wt.iter) {
				loc = "left"
			} else {
				loc = "topright"
			} 
		}
		legend(loc, bty = "n", paste(lvl, "%", c("Bootstrap C.I.  "), "[", round(ci[1], 4), ",", round(ci[2], 4), "]"), col = color, cex = 1)
		abline(v = ci[1], lty = 2)
		abline(v = ci[2], lty = 2)
	}
	if (has.leg) {
		if (is.null(zmat)) {
  			legend("topleft", bty = "n", c("f_hat"), col = color, lty = 2, lwd = 2)
		} else {
			legend(x = min(x0s), y = rmax , bty = "n", c("constrained f_hat"), col = color, lty = 2, lwd = 2)
			#if (!is_fac) {
			#	znm_3 = znm_2
			#} else {	
			#	znm_3 = paste("covariate the ", 1:lzb, "th column")
			#}
			#if (lzb >= 1) {
			#	for (i in 1:lzb) {	
			#		legend(x = min(x0s), y = rmax - i * .12 * rmax, bty = "n", znm_3[i], col = palette[i], lty = 1, lwd = 2)
			#	}
			#}
		}
  	}
  	#points(knots, 1:length(knots)*0 + rmin, pch = 'X')
  	#if (categ == "jp") {
        #	abline(v = sub[1], lty = 2, col = 2, lwd = 2)
    	#	abline(v = sub[2], lty = 2, col = 2, lwd = 2)
  	#}
}


###########################
#print.summary.ShapeChange#
###########################
print.summary.ShapeChange = function(x,...) {
	if (!is.null(x$zcoefs)) {
		cat("Call:\n")
		print(x$call)
		cat("\n")
		cat("Coefficients:")
		cat("\n")
		printCoefmat(x$coefficients, P.values = TRUE, has.Pvalue = TRUE)
		cat("\n")
	} else {
		print ("No predictor is defined")	
	}
}

#####################
#summary.ShapeChange#
#####################
#summary.ShapeChange = function(object,...) {
#	if (!is.null(object$zcoefs)) {
#		family = object$family 
#		wt.iter = object$wt.iter
#		coefs = object$zcoefs
#		se = object$se.beta
#		tval = object$tval
#    		pvalbeta = object$pvals.beta
#		n = length(coefs)
#		zmat = object$zmat
#		is_fac = object$is_fac
#		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)"
#		} 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:n] = colnames(zmat)
#		rslt1 = as.matrix(rslt1)
#		ans = list(call = object$call, coefficients = rslt1, zcoefs = coefs, family = family)
		#class(ans) = "summary.ShapeChange"
		#ans
#	} else {
#		ans = list(zcoefs = object$zcoefs)
		#class(ans) = "summary.ShapeChange"
		#ans
#	}
#	class(ans) = "summary.ShapeChange"
#	ans
#}

#####################
#predict.ShapeChange#
#####################
predict.ShapeChange = function(object, newData,...) {
	family = object$family
	cicfamily = CicFamily(family)
	muhat.fun = cicfamily$muhat.fun	
	coefs = object$coefs
	zcoefs = object$zcoefs
  	categ = object$categ
	sh = object$sh
  	tt = object$tms
 	knots = object$knots
  	nk = length(knots)
  	x = object$x
  	xs = sort(x)
  	xs = (xs - min(xs)) / (max(xs) - min(xs))
  	wt.iter = object$wt.iter
  	bmat = object$bmat
  	zmat = object$zmat
  	vals = object$vals
  	add = FALSE
	if (!inherits(object, "ShapeChange")) { 
    		warning("calling predict.ShapeChange(<fake-ShapeChange-object>) ...")
  	}
	if (missing(newData) | is.null(newData)) {
		ans = list(fhat = object$fhat, fhat_x = object$fhat_x)
		return (ans) 
	}
	Terms = delete.response(tt)
	m = model.frame(Terms, newData)
	newdata = m
  	nx = newdata[, 1]
  	if (!is.null(zcoefs)) {
     		nz = newdata[, 2]
     		if (!any(nz %in% vals)) {
        		stop ('New categorical variable value not found in the fit!')
     		}
#check!
     		if (length(unique(nz)) > 1) {
		        zmat = model.matrix(~ factor(nz))[, -1, drop = FALSE]
	     	} else {
	    		pos = which(vals == unique(nz))
       			add = TRUE
     		}
#print (head(zmat))
	}
  	nx0 = nx
  	snx = sort(nx)
#print (nx)
	nx = (snx - min(x)) / (max(x) - min(x))
  	if (any(nx0 > max(x)) | any(nx0 < min(x))) {
     		warning ('User is extrapolating!')
  	} else {
      		nx = (snx - min(x)) / (max(x) - min(x))
  	}
  	if (categ == 'mode' | categ == 'jp') {
     		ans0 = bqspl(x = nx, m = nk, knots = knots, pic = FALSE, x0 = xs)
  	} else if (categ == 'inflect') {
    		ans0 = bcspl(x = nx, m = nk, knots = knots, pic = FALSE, x0 = xs)
  	}
  	nbmat = ans0$bmat
	if (is.null(zcoefs)) {
     		nthb = nbmat %*% coefs
     		nthb_x = nthb
  	} else {
    		if (categ == 'mode' | categ == 'jp') {
       			nthb_x = nbmat[, 1:(nk + 1), drop = FALSE] %*% coefs[1:(nk + 1)]
    		} else if (categ == 'inflect') {
       			nthb_x = nbmat[, 1:(nk + 2), drop = FALSE] %*% coefs[1:(nk + 2)]
    		}
    		if (!add) {
       			nbmat = cbind(nbmat, zmat)
       			nthb = nbmat %*% coefs
    		} else {
      			if (pos == 1) {nthb = nthb_x} else {nthb = nthb_x + zcoefs[pos - 1]}
    		}
	}   
	if (wt.iter & categ == 'mode') {
     		if (is.null(zcoefs)) {
        		nthb = muhat.fun(nthb, fml = family$family)
        		nthb_x = nthb
     		} else {
       			nthb_x = muhat.fun(nthb_x, fml = family$family)
       			nthb = muhat.fun(nthb, fml = family$family)
     		}
  	}
	ans = list(fhat = nthb, fhat_x = nthb_x, nbmat = nbmat)
	return (ans) 
 }

############################
#find penalty terms #
############################
find_pen = function(aims = NULL, Q = NULL, B = NULL, D = NULL, PNT = TRUE, Y = NULL, D0 = NULL) {
        la = length(aims)
        lam_use = lambdas = NULL
	#gcvus = NULL
	if (PNT) {
            for (i in 1:la) {
              x.l = 1e+10; x.r = 1e-10                        
              #f = function(pen0) sum(diag(B %*% solve((Q + pen0 * D), t(B)))) - aims[i]
	      f = function(pen0) {
	      	umat = chol(Q + pen0 * D)
		#amat = B %*% solve(umat)
        	#umat = chol(qv)
        	iumat = diag(ncol(umat))
        	uinv = backsolve(umat, iumat)
		amat = B %*% uinv
		sum(amat * amat) - aims[i]
	      }
              if (abs(f(x.r)) < 1e-4) {
#print ('no pen')
                 lambda = 0
              } else {
                #lambda = uniroot(f, c(x.l, x.r), tol = 1e-8)$root
		lambda = uniroot(f, c(x.l, x.r))$root
              }
              lambdas = c(lambdas, lambda)
           }
        } else {
           lambdas = 0
	}
	return (lambdas)
}


#############
#mode.kts#
#############
mode.kts = function(x, y, zmat = NULL, m = NULL, knots = NULL, sh = 1, q = 3, pen = 0, pnt = FALSE, wt.iter = FALSE, arp = FALSE, lambdas = NULL, pen_bt = NULL, p_bt = NULL, hs = NULL, ids = NULL, family = gaussian(), gcv = FALSE, spl = NULL, dmat = NULL, x1 = NULL, xn = NULL) {
	linkfun = family$linkfun
	cicfamily = CicFamily(family)
	llh.fun = cicfamily$llh.fun
	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 
	xs = sort(x)	
	ord = order(x)
  y = y[ord]
#new:must order zmat too!
	if (!is.null(zmat)) {
		nzmat = apply(zmat, 2, function(elem) elem[ord])
		zmat = nzmat
	}
#if (!is.null(zmat)) {
#	nzmat = zmat
#	for (i in 1:ncol(zmat)) {
#		nzmat[, i] = (zmat[, i])[ord]
#
#	}
#	zmat = nzmat
#}
	if (sh == -1 & !wt.iter) {
		y = -y
	}	
	if (is.matrix(y)) {
		y = as.numeric(y)
	}
#new:z
	x = (xs - min(xs)) / (max(xs) - min(xs))
	n = length(x)	
	sm = 1e-5
	if (is.null(m)) {
		m0 = 3 + round(n^(1 / 7)) 
		m1 = m0
		m = NULL
		while (is.null(m)) { 
			pts = floor((n - 2) / (m1 - 1))
			rem_pts = (n - 2) %% (m1 - 1)
			if (pts > 2) {
				m = m1
			} else if (pts == 2) {
				if (rem_pts / (m1 - 1) >= 1) {
					m = m1
				} else {
					m1 = m1 - 1	
				}
			} else {
				m1 = m1 - 1
			}
		}
	}
	m2 = m
#new:
	if (is.null(spl)) {
		ans_spl = bqspl(x, m = m2, knots = knots, pnt = pnt)
	}		
	knots = ans_spl$knots
	bmat = ans_spl$bmat
  	slopes = ans_spl$slopes
#new:
	if (!is.null(zmat)) {
		bmat = cbind(bmat, zmat)
    		np = ncol(zmat)
		m = ncol(bmat) - np
	} else {
		np = 0
		m = ncol(bmat) 	
	}
  	qv0 = crossprod(bmat)
  	qv = qv0
  	#slopes = t(ans$slopes)
	la = 1
#new: return edfs
	dmat = edfs = gcvu = gcvus = lambda = lambdas = lambdas_pen = NULL
	if (is.null(dmat)) {
		#ans_dmat = make_dmat(m, q, np, pnt = pnt, wt.iter = wt.iter, arp = arp, gcv = gcv, p_bt = p_bt, type = 'mode')
		dmat = make_dmat(m, q, pnt)
	}
	#dmat = ans_dmat$dmat
#new: no pnt on zmat
if (pnt) {
	if (!is.null(zmat)) {
		dv0 = matrix(0, nrow = (m + np), ncol = (m + np))
		dv0[1:m, 1:m] = crossprod(dmat)
	} else {dv0 = crossprod(dmat)}
    	if (!wt.iter) {
       		if (!arp) {
			#if (is.null(pen)) {
			if (pen == 0L) {
				if (m < 6) {
					if (!gcv) {
# the unconstr number of columns of bmat 
						edfs = m
					} else {
						edfs = m:(m+2)
						if (!is.null(zmat)) {
							edfs = edfs + np
						}
					}
					#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
	           			pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					lambdas_pen = pen
				} else {
					if (!gcv) {
# the unconstr number of columns of bmat 
						edfs = 6
					} else {
						edfs = seq(4, m, by = 1)
						if (!is.null(zmat)) {
							edfs = edfs + np
						}
#print (edfs)
					}
					#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					lambdas_pen = pen 
				}
		        }
          		lambdas = pen
#print (lambdas)
       		}
       		if (arp & is.null(p_bt)) {        
	        	if (m < 10) {
            			edfs = 4:m
         		} else {
            			#aims = 5:10
         			edfs = seq(4, m, by = 1)
			}
			if (!is.null(zmat)) {
				edfs = edfs + np
			}
         		if (is.null(lambdas)) {
            			#if (is.null(pen)) {
				if (pen == 0L) {
					#ans_pen = find_pen(aims = aims, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
               				lambdas = pen
            			} else {
               				lambdas = pen
            			}
         		}
       		} 
      		if (!is.null(p_bt)) {
         		lambdas = as.vector(pen_bt)
       		}
       		la = length(lambdas)
    	}
}
#print (la)
	#lambdas = ans_dmat$lambdas
	#la = ans_dmat$la
	#dfmat = bmat %*% solve(qv, t(bmat))
  	pos = pos2 = NULL
	tr = trr = sig = aics = es = rmat = phi = minsse = NULL
	aiclst = aicmat = sigmat = trsmat = trrsmat = trslst = psmat = fsmat = phismat = psmat = pensmat = NULL
  	tr_use0 = tr_use = trr_use = NULL
#new:	
	edfc = edfcs = NULL
	id_gcv = gcvc = gcvcs = ps_gcv = ssegcv = ssegcvs = NULL
	fhats_gcv = fhats_x_gcv = bh_gcv = list()
	ms_gcv = 1:la*0
	tdf = ncol(bmat)
    	imat = diag(tdf)
	if (!wt.iter) {
		cv = crossprod(bmat, y)
    		smat = slopes
		if (!is.null(zmat)) {
			nilmat = matrix(0, nrow = nrow(smat), ncol = ncol(zmat))
			smat = cbind(smat, nilmat)
		}
#search for a window 
    		for (ila in 1:la) {
			#pen_ila = NULL
        		if (pnt) {
           			#pen_ila = lambdas[ila]
           			qv = qv0 + lambdas[ila] * dv0
        			#dfmat = bmat %*% solve(qv, t(bmat))
				#tr_use0 = c(tr_use0, sum(diag(dfmat)))
        		} else {qv = qv0}   
        		minsse = sum(y^2)
			#minsse = sum((y - mean(y))^2)
        		for (j in 1:m2) {
          			sl = smat 
          			sl[j:m2, ] = -sl[j:m2, ]    
				#bvec = 1:nrow(sl)*0      
          			qans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)
          			bh = fitted(qans)
#print (bh)
#new: get constrained edf
				qdf = qans$df
				d_use = t(qans$xmat)
				if (qdf < tdf) {
		             		pd = d_use %*% solve(crossprod(d_use), t(d_use))
		             		pmat0 = imat - pd
		          	} else {
		            		pmat0 = imat
		          	}             
				umat = chol(qv)
				#uinv = solve(umat)
          			iumat = diag(ncol(umat))
          			uinv = backsolve(umat, iumat)
          			#bu = bmat %*% uinv    
          			#pmat = bu %*% tcrossprod(pmat0, bu)
				pmat = bmat %*% uinv %*% pmat0 %*% t(uinv) %*% t(bmat)
          			fhat = bmat %*% bh
		         	edfcj = sum(diag(pmat))
#print (c(j, edfcj))
				ssegcvj = sum((y - fhat)^2) 
				gcvcj = ssegcvj / (1 - edfcj/n)^2
          			fhat_x = bmat[, 1:(m2 + 1), drop = FALSE] %*% bh[1:(m2 + 1), , drop = FALSE]
				sse = mean((y - fhat)^2)
#print (c(j, gcvcj))
				#sse = sum((y - fhat)^2) 
#print (sse)
				#if (sse < minsse) {thb = fhat; thb_x = fhat_x; minsse = sse; pos = j; bhat = bh; tr = tr_use0; lambda = lambdas[ila]; edfc = edfcj; gcvc = gcvcj; ssegcv = ssegcvj}
				if (sse < minsse) {thb = fhat; thb_x = fhat_x; minsse = sse; pos = j; bhat = bh; edfc = edfcj; gcvc = gcvcj; ssegcv = ssegcvj}
#print (edfc)
       		      	} 
			pos2 = c(pos2, pos)
       			pos_use = round(mean(pos2))
#new:
			edfcs = c(edfcs, edfc)
#print (edfcs)
			gcvcs = c(gcvcs, gcvc)
			ssegcvs = c(ssegcvs, ssegcv)
			fhats_gcv[[ila]] = thb
			fhats_x_gcv[[ila]] = thb_x
			bh_gcv[[ila]] = bhat
			ps_gcv = c(ps_gcv, pos)
		} 
    		if (arp) {
        		st = pos_use - 2
        		ed = pos_use + 2   
			#st = 1
			#ed = m2        		
			if (st < 1) st = 1
        		if (ed > m2) ed = m2
			#minllh = 1e+100 #sum(y^2)
			minaic = sum(y^2)			
			aiclst = fslst = siglst = phislst = penslst = trslst = trrslst = list()
			iterlst = 0
			miniter = NULL
        		for (j in st:ed) {
				iterlst = iterlst + 1
            			sl = smat 
		      		sl[j:m2, ] = -sl[j:m2, ]
            			qans = make_arp(y, pnt = pnt, dmat = dmat, bmat = bmat, sl = sl, p_max = 2, m2 = m2, lambdas = lambdas, pen_fit = pen_bt, p_fit = p_bt, dv0 = dv0, qv0 = qv0, cv = cv, hs = hs, ids = ids)
				bh = qans$thetahat
            			theta = qans$fhat_use
				theta_x = bmat[, 1:(m2 + 1), drop = FALSE] %*% bh[1:(m2 + 1), , drop = FALSE] 
				aiclst[[iterlst]] = qans$aics	
				siglst[[iterlst]] = qans$sigs
				trslst[[iterlst]] = qans$trs
				trrslst[[iterlst]] = qans$trrs
				fslst[[iterlst]] = qans$fhats
				phislst[[iterlst]] = qans$phis
				penslst[[iterlst]] = qans$pens 
            			r_use = qans$r_use
				phi_use = qans$phi_use
            			sig_use = qans$sig_use
            			aics_use = qans$aics
            			tr_use = qans$tr_use
				trr_use = qans$trr_use
            			es_use = qans$es_use
            			lambda_use = qans$lambda_use       								
				aicj = qans$aicmin
	          		if (aicj < minaic) {thb = theta; thb_x = theta_x; minaic = aicj; pos = j; bhat = bh; miniter = iterlst; phi = phi_use; sig = sig_use; aics = aics_use; tr = tr_use; trr = trr_use; es = es_use; lambda = lambda_use; rmat = r_use}
        		}         
			aicmat = aiclst[[miniter]]
			sigmat = siglst[[miniter]] 
			trsmat = trslst[[miniter]]
			trrsmat = trrslst[[miniter]]  
			fsmat = fslst[[miniter]]   
			phismat = phislst[[miniter]] 
			pensmat = penslst[[miniter]]  
			id_minaic = which(aicmat == min(aicmat))
			aics = aicmat 			                         		
		}
	} else {
		if (!is.null(zmat)) {
			nc = m2 + 1 + ncol(zmat)
			#np = ncol(zmat)
       		} else {nc = m2 + 1}
		smat = slopes
		if (!is.null(zmat)) {
			nilmat = matrix(0, nrow = nrow(smat), ncol = ncol(zmat))
			smat = cbind(smat, nilmat)
       		}
       		qv0 = crossprod(bmat)
       		qv = qv0
       		if (pnt) {       
#print (pnt)    
          		#dv0 = crossprod(dmat)
#new: no pnt on zmat
			if (!is.null(zmat)) {
				dv0 = matrix(0, nrow = (m + np), ncol = (m + np))
				dv0[1:m, 1:m] = crossprod(dmat)
			} else {dv0 = crossprod(dmat)}
          		#if (is.null(pen)) {
			if (pen == 0L) {
				#if (m < 6) {
				edfs0 = min(c(m, 6))	
					if (!gcv) {
# the unconstr number of columns of bmat 
						edfs = edfs0
					} else {
						edfs = edfs0:(edfs0+2)
						if (!is.null(zmat)) {
							edfs = edfs + np
						}
					}
					#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
             				pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					lambdas_pen = pen 
				#} else {
				#	if (!gcv) {
# the unconstr number of columns of bmat 
				#		edfs = 6
				#	} else {
				#		edfs = 6:8
				#		if (!is.null(zmat)) {
				#			edfs = edfs + np
				#		}
				#	}
				#	#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
				#	pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
				#	lambdas_pen = pen
				#}
	        	}
			lambda = pen; lambdas = pen 
          		qv = qv0 + pen * crossprod(dmat)
       		}
       		#dfmat = bmat %*% solve(qv, t(bmat))
       		#tr = sum(diag(dfmat))
		minllh = 1e+3
		for (j in 1:m2) {
			etahat = etahat.fun(n, y, fml = family$family)
			gr = gr.fun(y, etahat, weights = NULL, fml = family$family)  
			wt = wt.fun(etahat, n, weights = NULL, fml = family$family)
			cvec = wt * etahat - gr
			qvk = crossprod(bmat, diag(wt)) %*% bmat
           		if (pnt) {
				qvk = qvk + pen * dv0
           		}
			cveck = crossprod(bmat, cvec)			
			sl = smat 
			sl[j:m2, ] = -sl[j:m2, ]
			qans = qprog(qvk, cveck, sl, 1:m2*0, msg = FALSE)
			bh = qans$thetahat
			etahat = bmat %*% bh
			muhat = muhat.fun(etahat, fml = family$family)
			diff = 1
			nrep = 0
			while (diff > sm & nrep < 100) {
				oldmu = muhat	
				nrep = nrep + 1
				gr = gr.fun(y, etahat, weights = NULL, fml = family$family)	
				wt = wt.fun(etahat, n, weights = NULL, fml = family$family) 
				cvec = wt * etahat - gr
				qvk = crossprod(bmat, diag(wt)) %*% bmat
				if (pnt) {
					qvk = qvk + pen * dv0
#print (pen)
		            	}
				cveck = crossprod(bmat, cvec)
				qans = qprog(qvk, cveck, sl, 1:m2*0, msg = FALSE)
				bh = qans$thetahat
				etahat = bmat %*% bh
				muhat = muhat.fun(etahat, fml = family$family)
				diff = mean((muhat - oldmu)^2)	
         		}
			etahat_x = bmat[, 1:(m2 + 1), drop = FALSE] %*% bh[1:(m2 + 1), ,drop = FALSE]
			muhat_x = muhat.fun(etahat_x, fml = family$family)
			d_use = t(qans$xmat)
			qdf = qans$df
			tdf = ncol(bmat)
			imat = diag(tdf)
			if (qdf < tdf) {
   				pd = d_use %*% solve(crossprod(d_use), t(d_use))
   				pmat0 = imat - pd
			} else {
  				pmat0 = imat
			}   
			umat = chol(qvk)
			iumat = diag(ncol(umat))
			uinv = backsolve(umat, iumat)
			bu = bmat %*% uinv
			pmat = bu %*% tcrossprod(pmat0, bu) %*% diag(wt)
			tr = sum(diag(pmat))
         		llh = llh.fun(y, muhat, etahat, n, weights = NULL, fml = family$family)
			if (llh < minllh) {thb = muhat; thb_x = muhat_x; thb_eta = etahat; thb_eta_x = etahat_x; minllh = llh; pos = j; bhat = bh; tr = tr}
		}
	}
#new: gcv
  	dist = round(slopes %*% bhat[1:(m2 + 1)], 6)
	chpt = find_chpt(dist = dist, pos = pos, m2 = m2, knots = knots)
#print (chpt)
	if (sh == -1) {
		thb = -thb
		thb_x = -thb_x
		bhat = -bhat
	}
#wrong se's...
	if (!is.null(zmat)) {
		zcoefs = bhat[(m2 + 2):(m2 + 2 + np - 1)]
		df_obs = sum(abs(bhat) > 0)
		prior.w = 1:n*0 + 1
		w = diag(as.vector(prior.w / deriv.fun(thb, fml = family$family)))
		sse1 = sum(prior.w * (y - thb)^2)
     		if ((n - np - 1.5 * (df_obs - np)) <= 0) {
		 	  sdhat2 = sse1
     		} else {
			 sdhat2 = sse1 / (n - np - 1.5 * (df_obs - np))
     		}
     		if (wt.iter) {
			  se2 = solve(t(zmat) %*% w %*% zmat)
     		} else {
			  se2 = solve(t(zmat) %*% diag(prior.w) %*% zmat) * sdhat2
     		}		 			
     		se.beta = sqrt(diag(se2))
     		tstat = zcoefs / se.beta
     		pvals.beta = 2 * (1 - pt(abs(tstat),  n - np - 1.5 * (df_obs - np))) 
  	} else {zcoefs = pvals.beta = se.beta = tstat = NULL}
#if (gcv) {
#	ms_gcv = ms_gcv * (max(xs) - min(xs)) + min(xs)
#	id_gcv = which(gcvcs == min(gcvcs))
#	chpt = ms_gcv[id_gcv]
#	thb = fhats_gcv[[id_gcv]] 
#	thb_x = fhats_x_gcv[[id_gcv]] 
#	bhat = bh_gcv[[id_gcv]]
#	if (!is.null(zmat)) {
#		zcoefs = bhat[(m2 + 2):(m2 + 2 + np - 1)]
# 	} else {zcoefs = pvals.beta = se.beta = tstat = NULL}
#} else {
	chpt = chpt * (max(xs) - min(xs)) + min(xs)
#print (chpt)
#	ms_gcv = ms_gcv * (max(xs) - min(xs)) + min(xs)
#}
#new: 
#print (edfs)
#print (gcvcs)
edfu = edfc = NULL
if (gcv) {
	if (length(edfcs) > 1 & !is.null(edfcs)) {
	edfc = edfcs[which(gcvcs == min(gcvcs ))]
	edfu = edfs[which(gcvcs == min(gcvcs ))]
	}
}
#print (edfs)
#print (edfcs)
#print (edfc)
#print (edfu)
#print (lambda)
	knots = knots * (max(xs) - min(xs)) + min(xs)
	if (!wt.iter) {thb_eta_x = thb_x; thb_eta = thb}
	ans = list(fhat = thb, fhat_x = thb_x, fhat_eta = thb_eta, fhat_eta_x = thb_eta_x, bhat = bhat, sse = minsse, pos = pos, chpt = chpt, knots = knots, dmat = dmat, dist = dist, sub = sub, zcoefs = zcoefs, pvals.beta = pvals.beta, se.beta = se.beta, tval = tstat, bmat = bmat, phi = phi, sig = sig, aics = aics, aiclst = aiclst, es = es, lambda = lambda, lams = lambdas, edfs = edfs, gcvus = gcvus, gcvu = gcvu, rmat = rmat, aicmat = aicmat, sigmat = sigmat, trsmat = trsmat, trrsmat = trrsmat, psmat = psmat, fsmat = fsmat, phismat = phismat, pensmat = pensmat, edfcs = edfcs, edfc = edfc, edfu = edfu, gcvcs = gcvcs, ms_gcv = ms_gcv, fhats_gcv = fhats_gcv, fhats_x_gcv = fhats_x_gcv, ssegcvs = ssegcvs, id_gcv = id_gcv, spl = ans_spl, dmat = dmat, x1 = NULL, xn = NULL)
	ans 
}

########
#ip.kts#
########
ip.kts = function(x, y, zmat = NULL, m = NULL, knots = NULL, q = 3, pen = 0, pnt = FALSE, sh = 1, fir = FALSE, wt.iter = FALSE, arp = FALSE, lambdas = NULL, pen_bt = NULL, p_bt = NULL, hs = NULL, ids = NULL, family = gaussian(), gcv = FALSE, spl = NULL, dmat = NULL, x1 = NULL, xn = NULL) {
#print ('start')
	linkfun = family$linkfun
	cicfamily = CicFamily(family)
	llh.fun = cicfamily$llh.fun
	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  
	xs = sort(x)	
	ord = order(x)
	y = y[ord]
#new: must order zmat too!
	if (!is.null(zmat)) {
		#nzmat = zmat
		nzmat = apply(zmat, 2, function(elem) elem[ord])
		zmat = nzmat
	}
#new:
  	mult = max(xs) - min(xs)
	x = (xs - min(xs)) / (max(xs) - min(xs))
	n = length(x)
	#xu = unique(x)
	#sm = 1e-5
#print ('to get m')
	if (is.null(m)) {
		m0 = 4 + round(n^(1 / 9)) 
		m1 = m0
		m = NULL
		while (is.null(m)) { 
			pts = floor((n - 2) / (m1 - 1))
			rem_pts = (n - 2) %% (m1 - 1)
			if (pts > 2) {
				m = m1
			} else if (pts == 2) {
				if (rem_pts / (m1 - 1) >= 1) {
					m = m1
				} else {
					m1 = m1 - 1	
				}
			} else {
				m1 = m1 - 1
			}
		}
	}
	m2 = m
#print ('finish m')
#print ('to get spl')
	if (is.null(spl)) {
#print ('null')
		spl = bcspl(x, m = m2, knots = knots, pnt = pnt)
	}
	knots = spl$knots
	bmat = spl$bmat
	secder = spl$secder
#print ('finish')
	if (!is.null(zmat)) {
		bmat = cbind(bmat, zmat)
	    	np = ncol(zmat)
	} else {np = 0}
#new
	if (is.null(zmat)) {
		m = ncol(bmat)
	} else {m = ncol(bmat) - np}
# qv0 include zmat
	qv0 = crossprod(bmat)
  	qv = qv0 
  	la = 1	
	dmat = dv0 = edfs = gcvu = gcvus = lambda = lambdas = lambdas_pen = tru = tr_use0 = tr_use = trr_use = NULL 
	if (is.null(dmat)) {
		dmat = make_dmat(m, q, pnt = pnt)
	}
#print ('pass pnt!')
#tm = proc.time()
#	dfmat = bmat %*% solve(qv, t(bmat))
#	tru = tr_use = sum(diag(dfmat))
#umat = chol(qv)
#iumat = diag(ncol(umat))
#uinv = backsolve(umat, iumat)
#uinv = solve(umat)
#bu = bmat %*% uinv
#tru = sum(bu * bu)
#print (tru)
#print (proc.time() - tm)
#print ('pass dfmat!')
	#dfmat = bmat %*% solve(qv, t(bmat))
#new: no pnt on zmat
if (pnt) {
		if (!is.null(zmat)) {
		dv0 = matrix(0, nrow = (m + np), ncol = (m + np))
		dv0[1:m, 1:m] = crossprod(dmat)
		} else {dv0 = crossprod(dmat)}
    		if (!wt.iter) {
       			if (!arp) {
          			#if (pen == 0) {
             			#	pen = find_pen(aims = 8, Q = qv0, B = bmat, D = dv0, PNT = pnt)
          			#}
				#if (is.null(pen)) {
				if (pen == 0L) {
					if (m < 8) {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = m
						} else {
							edfs = m:(m+2)
						}
						#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
	             				pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen 
					} else {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = 8
						} else {
							edfs = 8:10
						}
						#ans_pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen 
					}
		        	}
          			lambdas = pen 
#print (lambdas)
#new: get gcvu anyway
				#pmatu = bmat %*% solve((qv0 + pen * dv0), t(bmat)) 
			 	#muhatu = pmatu %*% y
			 	#sseu = sum((y - muhatu)^2)
			 	#edfu = sum(diag(pmatu))
			 	#gcvu = sseu / (1 - edfu/n)^2
       			}
       			if (arp & is.null(p_bt)) {        
	       			if (m < 10) {
            				aims = 5:m
         			}  else {
            				aims = 5:10
         			}
				edfs = aims
         			if (is.null(lambdas)) {
					pen = find_pen(aims = aims, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
            				lambdas = pen
         			}
       			} 
       			if (!is.null(p_bt)) {
          			lambdas = as.vector(pen_bt)
       			}
       			la = length(lambdas)
    	}
}
	aics = minsse = nrep = pos = pos2 = tr = sig = aics = es = rmat = phi = NULL
	aicmat = sigmat = trsmat = trrsmat = trslst = psmat = fsmat = phismat = psmat = pensmat = NULL
	if (!wt.iter) {
		if (sh == -1) {
#concave-convex
			y = -y
		}
		cv = crossprod(bmat, y)
		smat = secder 
		if (!is.null(zmat)) {
			nilmat = matrix(0, nrow = nrow(smat), ncol = ncol(zmat))
			smat = cbind(smat, nilmat)
		}
    		tdf = ncol(bmat)
    		imat = diag(tdf)
#new
#print ('to get 1st deri!')
		if (is.null(x1)) {
    			x1 = bcsplfirderi(x[1], knots)
		}
		if (is.null(xn)) {
    			xn = bcsplfirderi(x[n], knots) 
		}
    		for (ila in 1:la) {
        		if (pnt) {
           			pen = lambdas[ila]
           			qv = qv0 + pen * dv0
				#dfmat = bmat %*% solve(qv, t(bmat))
				#tr_use0 = c(tr_use0, sum(diag(dfmat)))
        		} else {qv = qv0}   
        		minsse = sum(y^2)
#ptm = proc.time()
		    	for (j in 1:m2) {
#print (j)
            			sl = smat 
				sl[j:m2, ] = -sl[j:m2, ]
			      	if (fir) {
#increase if sh == 1; decrease if sh == -1 
					row_add = matrix(1:(2 * (ncol(sl))) * 0, nrow = 2)
				       	sl = rbind(sl, row_add)
               				sl[m2 + 1, 1:(m2 + 2)] = x1 
               				sl[m2 + 2, 1:(m2 + 2)] = xn 
            			}
			      	qans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)
            			bh = fitted(qans)
            			qdf = qans$df
            			fhat = bmat %*% bh
            			fhat_x = bmat[, 1:(m2 + 2), drop = FALSE] %*% bh[1:(m2 + 2), ,drop = FALSE]
#lines(x, fhat_x, col = ila, lty=2)
            			sse = sum((y - fhat)^2)
			      	if (sse < minsse) {thb = fhat; thb_x = fhat_x; minsse = sse; pos = j; bhat = bh; tr = tr_use0; lambda = pen}   
       			}
       			pos2 = c(pos2, pos)
       			pos_use = round(mean(pos2))
    		}
#print ('loop')
#print (proc.time() - ptm)
   	 	if (arp) {
        		st = pos_use - 2
        		ed = pos_use + 2   
        		if (st < 1) st = 1
        		if (ed > m2) ed = m2
        		minllh = sum(y^2) #1e+3
			aiclst = fslst = siglst = pslst = phislst = penslst = trslst = trrslst = list()
			iterlst = 0
			miniter = NULL
        		for (j in st:ed) {
				iterlst = iterlst + 1
            			sl = smat 
			      	sl[j:m2, ] = -sl[j:m2, ]
            			if (fir) {
#increase if sh == 1; decrease if sh == -1 
               				row_add = matrix(1:(2 * (ncol(sl))) * 0, nrow = 2)
				       	sl = rbind(sl, row_add)
               				sl[m2 + 1, 1:(m2 + 2)] = x1 
               				sl[m2 + 2, 1:(m2 + 2)] = xn  
           			}
           			qans = make_arp(y, pnt = pnt, dmat = dmat, bmat = bmat, sl = sl, p_max = 2, m2 = m2, lambdas = lambdas, pen_fit = pen_bt, p_fit = p_bt, dv0 = dv0, qv0 = qv0, cv = cv, hs = hs, ids = ids)
				bh = qans$thetahat
            			theta = qans$fhat_use
				theta_x = bmat[, 1:(m2 + 2), drop = FALSE] %*% bh[1:(m2 + 2), , drop = FALSE] 
				aiclst[[iterlst]] = qans$aics	
				siglst[[iterlst]] = qans$sigs
				trslst[[iterlst]] = qans$trs
				trrslst[[iterlst]] = qans$trrs
				fslst[[iterlst]] = qans$fhats
				phislst[[iterlst]] = qans$phis
				penslst[[iterlst]] = qans$pens 
            			r_use = qans$r_use
				phi_use = qans$phi_use
            			sig_use = qans$sig_use
            			aics_use = qans$aics
            			tr_use = qans$tr_use
				trr_use = qans$trr_use
            			es_use = qans$es_use
            			lambda_use = qans$lambda_use          								
				llh = crossprod((y - theta), chol2inv(chol(r_use))) %*% (y - theta)					
				if (llh < minllh) {thb = theta; thb_x = theta_x; minllh = llh; pos = j; bhat = bh; miniter = iterlst; phi = phi_use; sig = sig_use; aics = aics_use; tr = tr_use; trr = trr_use; es = es_use; lambda = lambda_use; rmat = r_use}
        		}                                
#print (ar(es, order.max=2))
			aicmat = aiclst[[miniter]]
			sigmat = siglst[[miniter]] 
			trsmat = trslst[[miniter]]
			trrsmat = trrslst[[miniter]] 
			fsmat = fslst[[miniter]]   
			phismat = phislst[[miniter]] 
			pensmat = penslst[[miniter]]  
			id_minaic = which(aicmat == min(aicmat))
			aics = aicmat	                 
    		}
	} else {
#concave-convex
        	if (pnt) {
			if (pen == 0L) {
				if (m < 8) {
					if (!gcv) {
# the unconstr number of columns of bmat 
						edfs = m
					} else {
						edfs = m:(m+2)
					}
					pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
        				#pen = ans_pen$lam_use
					#gcvus = ans_pen$gcvus
					lambdas_pen = pen
				} else {
					if (!gcv) {
# the unconstr number of columns of bmat 
						edfs = 8
					} else {
						edfs = 8:10
					}
					pen = find_pen(aims = edfs, Q = qv0, B = bmat, D = dv0, PNT = pnt, Y = y, D0 = dmat)
					#pen = ans_pen$lam_use
					#gcvus = ans_pen$gcvus
					lambdas_pen = pen
				}
			} 
			qv = qv0 + pen * dv0
			lambda = lambdas = pen
#new: get gcvu anywa
			pmatu = bmat %*% solve((qv0 + pen * dv0), t(bmat)) 
			muhatu = pmatu %*% y
			sseu = sum((y - muhatu)^2)
			edfu = sum(diag(pmatu))
			gcvu = sseu / (1 - edfu/n)^2
         	} 
         	#dfmat = bmat %*% solve(qv, t(bmat))
         	#tr_use = sum(diag(dfmat))
         	#tru = tr_use
                tdf = ncol(bmat)
                imat = diag(tdf)
         	if (family$family == "binomial") {
# use solve.QP if the last b > 1
# no zmat; no vmat
              		#cv = crossprod(bmat, y)
		        smat = secder 
              		#minsse = sum(y^2)
              		minllh = 1e+3
		        for (j in 1:m2) {
#print (j)
				sl = smat 
			        sl[j:m2, ] = -sl[j:m2, ]            
                  		#if (sh == 1) {
                     		row_add = matrix(1:(3 * ncol(sl)) * 0, nrow = 3)
                     		sl = rbind(sl, row_add)
                     		sl[m2 + 1, 1] = 1
				if (is.null(x1)) {
		    			x1 = bcsplfirderi(x[1], knots)
				}
				if (is.null(xn)) {
		    			xn = bcsplfirderi(x[n], knots) 
				}
                     		#sl[m2 + 2, 1:(m2 + 2)] = bcsplfirderi(x[1], knots = knots)
                     		#sl[m2 + 3, 1:(m2 + 2)] = bcsplfirderi(x[n], knots = knots)
                     		sl[m2 + 2, 1:(m2 + 2)] = x1
                     		sl[m2 + 3, 1:(m2 + 2)] = xn
                  		#}
                  		diff = 1
                  		nrep = 0 
                  		sm = 1e-5
                  		muhat = 1:n*0 + .5
                  		wt = 1 / muhat / (1 - muhat)
                  		while (diff > sm & nrep < 100){
#print (nrep)
                        		oldmu = muhat	
                        		nrep = nrep + 1										
                        		if (pnt) {
                           			qv = crossprod(bmat, diag(wt)) %*% bmat + pen * dv0    
                        		} else {
                          			qv = crossprod(bmat, diag(wt)) %*% bmat
                        		}
                        		cv = crossprod(bmat, diag(wt)) %*% y
                        		ans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)
                        		bh = fitted(ans)
					if (bh[m2 + 2] > 1) {
                           			B = matrix(1:ncol(sl)*0, nrow = 1, ncol = ncol(sl))  
                           			B[m2 + 2] = 1
                           			b0 = qr.solve(qr(B), 1)
                           			H = qr.Q(qr(t(B)), complete = TRUE)[, -(1:(qr(t(B))$rank)), drop = FALSE]	
                           			atil = sl %*% H
                           			bvec = -sl %*% b0
                           			if (pnt) {
							qv1 = t(H) %*% (crossprod(bmat, diag(wt)) %*% bmat + pen * dv0) %*% H 
                           			} else {
                             				qv1 = t(H) %*% crossprod(bmat, diag(wt)) %*% bmat %*% H
                           			}
                           			cv1 = t(H) %*% crossprod(bmat, diag(wt)) %*% (y - bmat %*% b0)
#new: scale qv1 and cv1				
						const = norm(qv1, "2")
                           			qans = solve.QP(qv1 / const, cv1 / const, t(atil), bvec)
                           			phi = qans$solution
                           			btil = H %*% phi                        			
						bh = btil + b0
                        		}   
                        		muhat = bmat %*% bh
	                      		diff = mean((muhat - oldmu)^2)	
		                    	if (any(round(muhat, 8) == 0)| any(round(muhat, 8) == 1)) {
			                     wt = rep(1e+2, n)
			                     id = round(muhat, 8) > 0 & round(muhat, 8) < 1
			                     wt[id] = 1 / muhat[id] / (1 - muhat[id])
                        		} else {wt = as.vector(1 / muhat / (1 - muhat))} 
                  		}
#tr when mu_n > 1???
                     		d_use = t(ans$xmat)
                     		qdf = ans$df
                     		if (qdf < tdf) {
                       		 	pd = d_use %*% solve(crossprod(d_use), t(d_use))
                        		pmat0 = imat - pd
                     		} else {
                       			pmat0 = imat
                     		} 
                     		umat = chol(qv)
                     		iumat = diag(ncol(umat))
                     		uinv = backsolve(umat, iumat)
                     		bu = bmat %*% uinv
                     		pmat = bu %*% tcrossprod(pmat0, bu) %*% diag(wt)   
                     		tr_use = sum(diag(pmat))  
                  		muhat[round(muhat, 8) == 0] = sm
		        	muhat[round(muhat, 8) == 1] = 1
                  		llh = llh.fun(y, muhat, etahat = NULL, n, weights = NULL, fml = family$family)
                  		if (llh < minllh) {thb = muhat; thb_x = muhat; minllh = llh; pos = j; bhat = bh; tr = tr_use}
             		}               
      		}
      		if (family$family == "poisson") {
         		smat = secder 
	 		minllh = 1e+3
         		for (j in 1:m2) {
				if (is.null(x1)) {
		    			x1 = bcsplfirderi(x[1], knots)
				}
				if (is.null(xn)) {
		    			xn = bcsplfirderi(x[n], knots) 
				}
             			sl = smat 
             			sl[j:m2, ] = -sl[j:m2, ] 
             			#if (sh == 1) {
             			row_add = matrix(1:(3 * ncol(sl)) * 0, nrow = 3)
             			sl = rbind(sl, row_add)
             			sl[m2 + 1, 1] = 1
             			#sl[m2 + 2, 1:(m2 + 2)] = bcsplfirderi(x[1], knots = knots)
				sl[m2 + 2, 1:(m2 + 2)] = x1
             			if (fir) {
             			#	sl[m2 + 3, 1:(m2 + 2)] = bcsplfirderi(x[n], knots = knots)
					sl[m2 + 3, 1:(m2 + 2)] = xn
             			} else {
             				sl[m2 + 3, m2 + 2] = 1
            			} 
             			#}
             			diff = 1
             			nrep = 0 
             			sm = 1e-5
             			muhat = 1:n*0 + mean(y)
             			etahat = 1:n*0
             			wt = 1 / muhat 
             			while (diff > sm & nrep < 100){
                   			oldmu = muhat	
                   			nrep = nrep + 1										
                  			 if (pnt) {
                      				qv = crossprod(bmat, diag(wt)) %*% bmat + pen * dv0    
                   			} else {
                     				qv = crossprod(bmat, diag(wt)) %*% bmat
                   			}
                   			cv = crossprod(bmat, diag(wt)) %*% y
                   			ans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)
                   			bh = fitted(ans)
                   			muhat = bmat %*% bh
                   			diff = mean((muhat - oldmu)^2)	
                   			if (any(round(muhat, 8) == 0)) {
                     				wt = rep(1e+2, n)
                      				id = round(muhat, 8) > 0
                      				wt[id] = 1 / muhat[id]
                   			} else {wt = as.vector(1 / muhat)} 
            			 }
             			d_use = t(ans$xmat)
             			qdf = ans$df
            			if (qdf < tdf) {
               				pd = d_use %*% solve(crossprod(d_use), t(d_use))
               				pmat0 = imat - pd
             			} else {
               				pmat0 = imat
             			} 
             			umat = chol(qv)
             			iumat = diag(ncol(umat))
             			uinv = backsolve(umat, iumat)
             			bu = bmat %*% uinv
             			pmat = bu %*% tcrossprod(pmat0, bu) %*% diag(wt)   
             			tr_use = sum(diag(pmat))  
             			muhat[round(muhat, 8) == 0] = sm
             			etahat = log(muhat)
             			llh = llh.fun(y, muhat, etahat, n, weights = NULL, fml = family$family)
             			if (llh < minllh) {thb = muhat; thb_x = muhat; minllh = llh; pos = j; bhat = bh; tr = tr_use}
        		}   
        		thb[round(thb, 8) == 0] = 0                    
      		}
	}
#print ('To get chpt!')
#ptm = proc.time()
  	dist = round(secder %*% bhat[1:(m2 + 2), ,drop = FALSE], 6) 
	chpt = find_chpt(dist = dist, pos = pos, m2 = m2, knots = knots)
   	if (sh == -1 & !wt.iter) {
		thb = -thb
		thb_x = -thb_x
		bhat = -bhat
   	}
#print (proc.time() - ptm) 
#wrong:
	if (!is.null(zmat)) {
		zcoefs = bhat[(m2 + 3):(m2 + 3 + np - 1)]
     		df_obs = sum(abs(bhat) > 0)
		#df_obs = df		
		prior.w = 1:n*0 + 1
 		w = diag(as.vector(prior.w / deriv.fun(thb, fml = family$family)))
     		sse1 = sum(prior.w * (y - thb)^2)
     		if ((n - np - 1.5 * (df_obs - np)) <= 0) {
			sdhat2 = sse1
     		} else {
			sdhat2 = sse1 / (n - np - 1.5 * (df_obs - np))
     		}
     		if (wt.iter) {
			se2 = solve(t(zmat) %*% w %*% zmat)
     		} else {
			se2 = solve(t(zmat) %*% diag(prior.w) %*% zmat) * sdhat2
     		}		 			
     		se.beta = sqrt(diag(se2))
     		tstat = zcoefs / se.beta
     		pvals.beta = 2 * (1 - pt(abs(tstat),  n - np - 1.5 * (df_obs - np))) 
  	} else {zcoefs = pvals.beta = se.beta = tstat = NULL}
	chpt = chpt * (max(xs) - min(xs)) + min(xs)
#new: scale knots back
	knots = knots * (max(xs) - min(xs)) + min(xs)
#print (proc.time() - ptm)
	ans = list(fhat = thb, fhat_x = thb_x, bhat = bhat, sse = minsse, chpt = chpt, knots = knots, pos = pos, dist = dist, sub = sub, zcoefs = zcoefs, pvals.beta = pvals.beta, se.beta = se.beta, tval = tstat, bmat = bmat, phi = phi, sig = sig, aics = aics, es = es, lambda = lambda, lams = lambdas, edfs = edfs, gcvus = gcvus, gcvu = gcvu, lambdas_pen = lambdas_pen, rmat = rmat, aicmat = aicmat, sigmat = sigmat, trsmat = trsmat, trrsmat = trrsmat, fsmat = fsmat, phismat = phismat, pensmat = pensmat, spl = spl, dmat = dmat, x1 = x1, xn = xn)
	 ans 
}

###########
#jp.pts   #
###########		
jp.pts = function(x, y, zmat = NULL, jpt, i = NULL, m = NULL, knots = NULL, q = 3, pen = 0, pnt = FALSE, up = TRUE, trd1 = -1, trd2 = -1, constr = TRUE, wt.iter = FALSE, arp = FALSE, lambdas = NULL, pen_bt = NULL, p_bt = NULL, hs = NULL, ids = NULL, family = gaussian(), gcv = FALSE) {	
	linkfun = family$linkfun
	cicfamily = CicFamily(family)
	llh.fun = cicfamily$llh.fun
	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 	
	sm = 1e-5
	if (!up) {
		y = -y; trd1 = -trd1; trd2 = -trd2
	}
	xs = sort(x)	
	ord = order(x)
	y = y[ord]	
	n = length(x)
#new:must order zmat too!
if (!is.null(zmat)) {
	nzmat = zmat
	for (i in 1:ncol(zmat)) {
		nzmat[, i] = (zmat[, i])[ord]

	}
	zmat = nzmat
}
#new:
  	mult = (max(xs) - min(xs))
	x = (xs - min(xs)) / (max(xs) - min(xs))
	n = length(x)
	#xu = unique(x)
	if (is.null(m)) {
		m0 = 4 + round(n^(1 / 9)) 
		m1 = m0
		m = NULL
		while (is.null(m)) { 
			pts = floor((n - 2) / (m1 - 1))
			rem_pts = (n - 2) %% (m1 - 1)
			if (pts > 2) {
				m = m1
			} else if (pts == 2) {
				if (rem_pts / (m1 - 1) >= 1) {
					m = m1
				} else {
					m1 = m1 - 1	
				}
			} else {
				m1 = m1 - 1
			}
		}
	}
	ans = bqspl(x, m = m, knots = knots, pic = FALSE, pnt = pnt)
	knots = ans$knots
	slopes = ans$slopes
	bmat = ans$bmat
	m = length(knots)
	kobs = 1:m
	m0 = length(bmat) / n 
	mj = m0 + 2
	obs = 1:n
	#jp is x[i]
	i = obs[xs == jpt]	
	#i = obs[round(xs, 10) == jpt]
#print (i)
#if jp is at some x[i], then there will be tiny increase
	#jp = (x[i]+x[i+1])/2
	#jp0 = x[i]
	djump = matrix(0, nrow = n, ncol = mj)
	djump[, 1:(mj - 2)] = bmat
	djump[1:i, mj - 1] = (2 * i + 1) / 2 - n; djump[(i + 1):n, mj - 1] = (2 * i + 1) / 2
	djump[1:i, mj] = 0; djump[(i + 1):n, mj] = x[x > (x[i] + x[i + 1]) / 2] - (x[i] + x[i + 1] ) / 2
	m_i = mean(djump[, mj])
	djump[, mj] = djump[, mj] - mean(djump[, mj])
#new:
	if (!is.null(zmat)) {
		djump = cbind(djump, zmat)
		np = ncol(zmat)
	} else {np = 0}
#new: pnt, qv
	gcvu = gcvus = NULL
	edfs = NULL
	lambda = lambdas = lambdas_pen = NULL
	tr = tru = sub = kts_in = NULL
	rmat = phi = sig = aics = tr = trr = es = theta = etahat = theta_x = etahat_x = NULL
	if (constr) {
		dist = knots - (x[i] + x[i + 1]) / 2
		obs = 1:m
		pos = NULL
		if (any(round(dist, 8) == 0)) {
			pos = min(obs[round(dist, 8) == 0])
			sub = c(knots[pos-1], knots[pos+1])
		} else {
			for (im in 1:(m - 1)) {
				if (dist[im] < 0 & dist[im+1] > 0) {
					pos = im
					break 
				}
			}
			sub = c(knots[pos], knots[pos+1])
		}
		bool = knots >= sub[1] & knots <= sub[2]
		kts_in = knots[bool]
		id_in = kobs[bool]
		kn = 1:(m + 3) < 0
		bool1 = knots > (x[i]+x[i+1])/2 & knots <= sub[2]		
		kn[1:m] = bool1

		smat = matrix(0, nrow = m + 3, ncol = mj)
		if (trd1 == -1 & trd2 == -1) {
			smat[id_in, 1:(mj - 2)] = -slopes[id_in, ]
#useless case:
		} else if (trd1 == 1 & trd2 == 1) {
			smat[id_in, 1:(mj - 2)] = slopes[id_in, ]
		} else if (trd1 * trd2 < 0) {
			dist = knots - (x[i] + x[i+1]) / 2
			obs = 1:m
			pos = NULL
			if (any(round(dist, 8) == 0)) {
				pos = min(obs[round(dist, 8) == 0])
			} else {
				for (im in 1:(m - 1)) {
					if (dist[im] < 0 & dist[im+1] > 0) {
						pos = im
						break 
					}
				}
			}
			bool2 = knots >= sub[1] & knots <= knots[pos]
			kts_in2 = knots[bool2]
			id_in2 = kobs[bool2]

			bool3 = knots > knots[pos] & knots <= sub[2]
			kts_in3 = knots[bool3]
			id_in3 = kobs[bool3]

			smat[id_in, 1:(mj - 2)] = slopes[id_in, ]
			if (trd1 == - 1 & trd2 == 1) {		
				smat[id_in2, 1:(mj - 2)] = -smat[id_in2, 1:(mj - 2)]		
			} 
			if (trd1 == 1 & trd2 == -1) {
				smat[id_in3, 1:(mj - 2)] = -smat[id_in3, 1:(mj - 2)]		
			}
		}
smat[m + 1, mj - 1] = 1
		if (trd1 == -1 & trd2 == -1) {
			smat[, mj] = 0
			if (any(bool1)) {
				smat[kn, mj] = -1
			} 
			smat[m + 2, mj] = -1
			ansi = sl_fun((x[i] + x[i+1]) / 2, knots, slopes)
			smat[m + 2, 1:m0] = -ansi 
			smat[m + 3, 1:m0] = -ansi  
		} else if (trd1 == 1 & trd2 == 1) {
			smat[, mj] = 0
			if (any(bool1)) {
				smat[kn, mj] = 1
			} 
			smat[m + 2, mj] = 1
			ansi = sl_fun((x[i] + x[i+1]) / 2, knots, slopes)
			smat[m + 2, 1:m0] = ansi 
			smat[m + 3, 1:m0] = ansi 
		} else if (trd1 == -1 & trd2 == 1) {
			smat[, mj] = 0
			if (any(bool1)) {
				smat[kn, mj] = 1
			} 
			smat[m + 2, mj] = 1
			ansi = sl_fun((x[i] + x[i+1]) / 2, knots, slopes)
			smat[m + 2, 1:m0] = ansi 
			smat[m + 3, 1:m0] = -ansi 
		} else if (trd1 == 1 & trd2 == -1) {
#smat[m + 1, mj - 1] = 1
			smat[, mj] = 0
			if (any(bool1)) {
				smat[kn, mj] = -1
			} 
			smat[m + 2, mj] = -1
			ansi = sl_fun((x[i] + x[i+1]) / 2, knots, slopes)
			smat[m + 2, 1:m0] = -ansi 
			smat[m + 3, 1:m0] = ansi 
		} 
		use = 1:(m + 3) > 0
		kp = min(knots[knots > (x[i]+x[i+1])/2])
		if (sum(x > (x[i]+x[i+1])/2 & x <= kp) == 0){use[m + 2] = FALSE}
		kp = max(knots[knots < (x[i]+x[i+1])/2])
		if (sum(x < (x[i]+x[i+1])/2 & x >= kp) == 0){use[m + 3] = FALSE}
		smat = smat[use, ]
		if (i == 1 | i == (n - 1)) {
			smat = smat[, -mj, drop = FALSE]
			djump = djump[, -mj, drop = FALSE] 
		}	
		if (!is.null(zmat)) {
			nilmat = matrix(0, nrow = nrow(smat), ncol = ncol(zmat))
			smat = cbind(smat, nilmat)
			#np = ncol(zmat)
		} 
		qv0 = crossprod(djump)
		qv = qv0		
		if (pnt) {
			#mj2 = ncol(djump)
#new
if (is.null(zmat)) {
	mj2 = ncol(djump)
} else {mj2 = ncol(djump) - np}
			dmat = matrix(0, nrow = (mj2 - q), ncol = mj2)
#  third-order
			if (q == 3) {
				for (idm in 4:mj2) {
					dmat[idm - 3, idm - 3] = 1; dmat[idm - 3, idm - 2] = -3; dmat[idm - 3, idm - 1] = 3; dmat[idm - 3, idm] = -1
				}
			}
# second order
			if (q == 2) {
				for (idm in 3:mj2) {
					dmat[idm - 2, idm - 2] = 1; dmat[idm - 2, idm - 1] = -2; dmat[idm - 2, idm] = 1
				}
			}
# first order
			if (q == 1) {
				for (idm in 2:mj2) {
					dmat[idm - 1, idm - 1] = 1; dmat[idm - 1, idm] = -1
				}
			}
# zero order
			if (q == 0) {
				for (idm in 1:mj2) {
					dmat[idm, idm] = 1
				}
			}    
    			#dv0 = crossprod(dmat)
#new: no pnt on zmat
if (!is.null(zmat)) {
	dv0 = matrix(0, nrow = (mj2 + np), ncol = (mj2 + np))
	dv0[1:mj2, 1:mj2] = crossprod(dmat)
} else {dv0 = crossprod(dmat)}
       			if (!arp) {
          			#if (pen == 0) {
             			#	pen = find_pen(aims = 6, Q = qv0, B = djump, D = dv0, PNT = pnt)
          			#}
				#if (pen == 0) {
				#if (is.null(pen)) {
				if (pen == 0L) {
					if (mj2 < 8) {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = mj2
						} else {
							edfs = mj2:(mj2+2)
						}
						pen = find_pen(aims = edfs, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
	             				#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen
					} else {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = 8
						} else {
							edfs = 8:10
						}
						pen = find_pen(aims = edfs, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen
					}
		        	}
          			lambdas = pen
#new: get gcvu anyway
				pmatu = djump %*% solve((qv0 + pen * dv0), t(djump)) 
			 	muhatu = pmatu %*% y
			 	sseu = sum((y - muhatu)^2)
			 	edfu = sum(diag(pmatu))
			 	gcvu = sseu / (1 - edfu/n)^2

       			}
       			if (arp & is.null(p_bt)) {        
	        		if (mj2 < 10) {
					aims = 5:mj2
         			} else {
            				aims = 5:10
         			}
				edfs = aims
         			if (is.null(lambdas)) {
            				if (pen == 0) {
						pen = find_pen(aims = aims, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
               					lambdas = pen
            				} else {
               					lambdas = pen
            				}
         			}
       			} 
      			if (!is.null(p_bt)) {
          			lambdas = as.vector(pen_bt)
       			}
       			la = length(lambdas)
			qv = qv + pen * dv0
		} 
		if (!wt.iter) {
			cv = crossprod(djump, y)
			if (!arp) {
	         		dfmat = djump %*% solve(qv, t(djump))
	         		tr_use = sum(diag(dfmat))
	         		tru = tr_use
				fit = qprog(qv, cv, smat, 1:nrow(smat)*0, msg = FALSE)				
				bhat = fitted(fit)		
				theta = djump %*% bhat
				nd = ncol(djump)
				if (!is.null(zmat)) {
					nd = nd - ncol(zmat)
				}
				theta_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, ,drop = FALSE]
				df = fit$df
			} else {
				qans = make_arp(y, pnt = pnt, dmat = dmat, bmat = djump, sl = smat, p_max = 4, m2 = m, lambdas = lambdas, pen_fit = pen_bt, p_fit = p_bt, dv0 = dv0, qv0 = qv0, cv = cv, hs = hs, ids = ids)
				bhat = qans$thetahat
	            		theta = qans$fhat_use
				nd = ncol(djump)
				if (!is.null(zmat)) {
					nd = nd - ncol(zmat)
				}
				theta_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, , drop = FALSE]
	            		rmat = qans$r_use
				phi = qans$phi_use
	            		sig = qans$sig_use
#print (sig)
	            		aics = qans$aics
	            		tr = qans$tr_use
				trr = qans$trr_use
	            		es = qans$es_use
	            		lambda = qans$lambda_use
			}
		} else {
			etahat = etahat.fun(n, y, fml = family$family)
			gr = gr.fun(y, etahat, weights = NULL, fml = family$family)  
			wt = wt.fun(etahat, n, weights = NULL, fml = family$family)
			cvec = wt * etahat - gr
			qvk = crossprod(djump, diag(wt)) %*% djump
	        	if (pnt) {
				#if (pen == 0) {
	              		#	pen = find_pen(aims = 6, B = djump, D = dmat, PNT = pnt)
				#}
				#if (pen == 0) {
				#if (is.null(pen)) {
				if (pen == 0L) {
					if (mj2 < 8) {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = mj2
						} else {
							edfs = mj2:(mj2+2)
						}
						pen = find_pen(aims = edfs, Q = qvk, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
	             				#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen
					} else {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = 8
						} else {
							edfs = 8:10
						}
						pen = find_pen(aims = edfs, Q = qvk, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen

					}
		        	}
				qvk = qvk + pen * dv0
				lambda = pen; lambdas = pen        
#new: get gcvu anyway
				pmatu = djump %*% solve((qvk + pen * dv0), t(djump)) 
			 	muhatu = pmatu %*% y
			 	sseu = sum((y - muhatu)^2)
			 	edfu = sum(diag(pmatu))
			 	gcvu = sseu / (1 - edfu/n)^2
			}
			cveck = crossprod(djump, cvec)			
			qans = qprog(qvk, cveck, smat, 1:nrow(smat)*0, msg = FALSE)
			bhat = qans$thetahat
			etahat = djump %*% bhat
			muhat = muhat.fun(etahat, fml = family$family)
			diff = 1
			nrep = 0
			while (diff > sm & nrep < 100) {
				oldmu = muhat	
				nrep = nrep + 1
				gr = gr.fun(y, etahat, weights = NULL, fml = family$family)	
				wt = wt.fun(etahat, n, weights = NULL, fml = family$family) 
				cvec = wt * etahat - gr
				qvk = crossprod(djump, diag(wt)) %*% djump
				if (pnt) {
					qvk = qvk + pen * dv0
				}
				cveck = crossprod(djump, cvec)	
				qans = qprog(qvk, cveck, smat, 1:nrow(smat)*0, msg = FALSE)
				bhat = qans$thetahat
				etahat = djump %*% bhat
				nd = ncol(djump)
				if (!is.null(zmat)) {
					nd = nd - ncol(zmat)
				}
				muhat = muhat.fun(etahat, fml = family$family)
				diff = mean((muhat - oldmu)^2)	
#lines(x, muhat, col = nrep + 1, lty = 2)
	       		} 
	       		etahat_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, ,drop = FALSE]
	       		muhat_x = muhat.fun(etahat_x, fml = family$family)
			theta = muhat
			theta_x = muhat_x
		}
	} else {
		if (!wt.iter) {
			if (i == 1 | i == (n - 1)) {
				djump = djump[, -mj, drop = FALSE] 
			}	
			qv0 = crossprod(djump)
			qv = qv0
			cv = crossprod(djump, y)
			if (pnt) {
				#mj2 = ncol(djump)
#new
if (is.null(zmat)) {
	mj2 = ncol(djump)
} else {mj2 = ncol(djump) - np}
				dmat = matrix(0, nrow = (mj2 - q), ncol = mj2)
#  third-order
				if (q == 3) {
					for (idm in 4:mj2) {
						dmat[idm - 3, idm - 3] = 1; dmat[idm - 3, idm - 2] = -3; dmat[idm - 3, idm - 1] = 3; dmat[idm - 3, idm] = -1
					}
				}
# second order
				if (q == 2) {
					for (idm in 3:mj2) {
						dmat[idm - 2, idm - 2] = 1; dmat[idm - 2, idm - 1] = -2; dmat[idm - 2, idm] = 1
					}
				}
# first order
				if (q == 1) {
					for (idm in 2:mj2) {
						dmat[idm - 1, idm - 1] = 1; dmat[idm - 1, idm] = -1
					}
				}
# zero order
				if (q == 0) {
					for (idm in 1:mj2) {
						dmat[idm, idm] = 1
					}
				}    
    				#dv0 = crossprod(dmat)
#new: no pnt on zmat
if (!is.null(zmat)) {
	dv0 = matrix(0, nrow = (mj2 + np), ncol = (mj2 + np))
	dv0[1:mj2, 1:mj2] = crossprod(dmat)
} else {dv0 = crossprod(dmat)}
       				if (!arp) {
					#if (pen == 0) {
					#if (is.null(pen)) {
					if (pen == 0L) {
						if (mj2 < 8) {
							if (!gcv) {
# the unconstr number of columns of bmat 
								edfs = mj2
							} else {
								edfs = mj2:(mj2+2)
							}
							pen = find_pen(aims = edfs, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
	             					#pen = ans_pen$lam_use
							#gcvus = ans_pen$gcvus
							lambdas_pen = pen
						} else {
							if (!gcv) {
# the unconstr number of columns of bmat 
								edfs = 8
							} else {
								edfs = 8:10
							}
							pen = find_pen(aims = edfs, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
							#pen = ans_pen$lam_use
							#gcvus = ans_pen$gcvus
							lambdas_pen = pen
						}
		        		}
          				lambdas = pen
#new: get gcvu anyway
					pmatu = djump %*% solve((qv0 + pen * dv0), t(djump)) 
			 		muhatu = pmatu %*% y
			 		sseu = sum((y - muhatu)^2)
			 		edfu = sum(diag(pmatu))
			 		gcvu = sseu / (1 - edfu/n)^2
       				}
       				if (arp & is.null(p_bt)) {        
	        			if (mj2 < 10) {
						aims = 5:mj2
         				} else {
            					aims = 5:10
         				}
					edfs = aims
         				if (is.null(lambdas)) {
            					#if (pen == 0) {
						#if (is.null(pen)) {
						if (pen == 0L) {
							pen = find_pen(aims = aims, Q = qv0, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
               						lambdas = pen
            					} else {
               						lambdas = pen
            					}
         				}
       				} 
      				if (!is.null(p_bt)) {
          				lambdas = as.vector(pen_bt)
       				}
       				la = length(lambdas)
				qv = qv + pen * dv0
			} 	
			if (!arp) {
				bhat = solve(qv, cv)				
				#bhat = solve(crossprod(djump), crossprod(djump, y))
				theta = djump %*% bhat		
				#nd = ncol(djump)
				#if (!is.null(zmat)) {
				#	nd = nd - ncol(zmat)
				#}
				#theta_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, ,drop = FALSE]
				#df = n - ncol(djump)
			} else {
				qans = make_arp(y, pnt = pnt, dmat = dmat, bmat = djump, sl = NULL, p_max = 4, m2 = m, lambdas = lambdas, pen_fit = pen_bt, p_fit = p_bt, dv0 = dv0, qv0 = qv0, cv = cv, hs = hs, ids = ids, constr = FALSE)
				bhat = qans$thetahat
	            		theta = qans$fhat_use
				#nd = ncol(djump)
				#if (!is.null(zmat)) {
				#	nd = nd - ncol(zmat)
				#}
				#theta_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, , drop = FALSE]
	            		rmat = qans$r_use
				phi = qans$phi_use
	            		sig = qans$sig_use
	            		aics = qans$aics
	            		tr = qans$tr_use
				trr = qans$trr_use
	            		es = qans$es_use
	            		lambda = qans$lambda_use
			}
			nd = ncol(djump)
			if (!is.null(zmat)) {
				nd = nd - ncol(zmat)
			}
			theta_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, ,drop = FALSE]
			df = n - ncol(djump)
		} else {
			etahat = etahat.fun(n, y, fml = family$family)
			gr = gr.fun(y, etahat, weights = NULL, fml = family$family)  
			wt = wt.fun(etahat, n, weights = NULL, fml = family$family)
			cvec = wt * etahat - gr
			qvk = crossprod(djump, diag(wt)) %*% djump
        		if (pnt) {
				#if (pen == 0) {
       		 	      	#	pen = find_pen(aims = 6, B = djump, D = dmat, PNT = pnt)
				#}
				#if (pen == 0) {
				#if (is.null(pen)) {
				if (pen == 0L) {
					if (mj2 < 8) {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = mj2
						} else {
							edfs = mj2:(mj2+2)
						} 
						pen = find_pen(aims = edfs, Q = qvk, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
             					#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen
					} else {
						if (!gcv) {
# the unconstr number of columns of bmat 
							edfs = 8
						} else {
							edfs = 8:10
						}
						pen = find_pen(aims = edfs, Q = qvk, B = djump, D = dv0, PNT = pnt, Y = y, D0 = dmat)
						#pen = ans_pen$lam_use
						#gcvus = ans_pen$gcvus
						lambdas_pen = pen 
					}
	        		}
				qvk = qvk + pen * dv0
				lambda = pen; lambdas = pen
#new: get gcvu anyway
				pmatu = djump %*% solve((qvk + pen * dv0), t(djump)) 
			 	muhatu = pmatu %*% y
			 	sseu = sum((y - muhatu)^2)
			 	edfu = sum(diag(pmatu))
			 	gcvu = sseu / (1 - edfu/n)^2
       	 		}
			cveck = crossprod(djump, cvec)			
			bhat = solve(qvk, cveck)
			etahat = djump %*% bhat
			muhat = muhat.fun(etahat, fml = family$family)
			diff = 1
			nrep = 0
			while (diff > sm & nrep < 100) {
				oldmu = muhat	
				nrep = nrep + 1
				gr = gr.fun(y, etahat, weights = NULL, fml = family$family)	
				wt = wt.fun(etahat, n, weights = NULL, fml = family$family) 
				cvec = wt * etahat - gr
				qvk = crossprod(djump, diag(wt)) %*% djump
				if (pnt) {
					qvk = qvk + pen * dv0
				}
				cveck = crossprod(djump, cvec)	
				bhat = solve(qvk, cveck)
				etahat = djump %*% bhat
				nd = ncol(djump)
				if (!is.null(zmat)) {
					nd = nd - ncol(zmat)
				}
				muhat = muhat.fun(etahat, fml = family$family)
				diff = mean((muhat - oldmu)^2)	
       			} 
			etahat_x = djump[, 1:nd, drop = FALSE] %*% bhat[1:nd, ,drop = FALSE]
			muhat_x = muhat.fun(etahat_x, fml = family$family)
			theta = muhat
			theta_x = muhat_x
		}
	}
	if (!up) {
		theta = -theta
		theta_x = -theta_x
		bhat = -bhat
	}
#wrong: 
	if (!is.null(zmat)) {
		mj2 = ncol(djump)
		zcoefs = bhat[(mj2 - np + 1):mj2]
     		df_obs = sum(abs(bhat) > 0)
		#df_obs = df		
		prior.w = 1:n*0 + 1
 		w = diag(as.vector(prior.w / deriv.fun(theta, fml = family$family)))
     		sse1 = sum(prior.w * (y - theta)^2)
     		if ((n - np - 1.5 * (df_obs - np)) <= 0) {
			sdhat2 = sse1
     		} else {
			sdhat2 = sse1 / (n - np - 1.5 * (df_obs - np))
     		}
     		if (wt.iter) {
			se2 = solve(t(zmat) %*% w %*% zmat)
     		} else {
			se2 = solve(t(zmat) %*% diag(prior.w) %*% zmat) * sdhat2
     		}		 			
     		se.beta = sqrt(diag(se2))
     		tstat = zcoefs / se.beta
     		pvals.beta = 2 * (1 - pt(abs(tstat),  n - np - 1.5 * (df_obs - np))) 
  	} else {zcoefs = pvals.beta = se.beta = tstat = NULL}
	sse = sum((y - theta)^2)
	ans = list(pos = i, chpt = jpt, fhat = theta, fhat_x = theta_x, fhat_eta = etahat, fhat_eta_x = etahat_x, df = df, sse = sse, knots = knots, sub = sub, kts_in = kts_in, bhat = bhat, zcoefs = zcoefs, m_i = m_i, sub = sub, lams = lambdas, lambda = lambda, edfs = edfs, gcvus = gcvus, gcvu = gcvu, lambdas_pen = lambdas_pen, tr = tr, trr = trr, tru = tru, phi = phi, sig = sig, aics = aics, es = es, rmat = rmat)
	ans
}

############################################################
#make the first derivative at a jump point by interpolation#
############################################################
sl_fun = function(xinterp, knots, slopes) {
	nk = length(knots)
	obs = 1:nk
	dist = knots - xinterp
	id = NULL; id1 = NULL
	if (any(round(dist, 8) == 0)) {
		id = min(obs[round(dist, 8) == 0])
	} else {	
		for (i in 1:(nk - 1)) {
			if (dist[i] < 0 && dist[i + 1] > 0) {
				id = i
				break 
			}
		}
	}
	k1 = knots[id]
	k2 = knots[id + 1]
	sinterp = 1:(nk + 1)*0
	a = (xinterp - k1) / (k2 - k1) 
	yk1 = slopes[id, ]
	yk2 = slopes[id + 1, ]
	sinterp = a * (yk2 - yk1) + yk1
	sinterp
}

#######################
#bqspline basis matrix# 
#######################
bqspl = function(x, m = NULL, knots = NULL, pic = FALSE, spl = TRUE, x0 = NULL, pnt = FALSE) {
############
#make knots#  
############
	bmat = slopes = NULL
	xu = unique(x)
#knots are equal x quantile of length = m spaced on the support of min(x) to max(x)
  	if (!is.null(knots)) {
		m = length(knots)
	} else {
		if (!is.null(m)) { 
			#knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			knots = 0:(m-1)/(m-1)*(max(x)-min(x))+min(x)
			#if (pnt) {
			#	knots = 0:(m-1)/(m-1)*(max(xu)-min(xu))+min(xu)
			#} else {
			#	knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#}
		} else {
			m0 = 4 + round(n^(1 / 9)) 
			m1 = m0
			#m = NULL
			while (is.null(m)) { 
				pts = floor((n - 2) / (m1 - 1))
				rem_pts = (n - 2) %% (m1 - 1)
				if (pts > 3) {
					m = m1
				} else if (pts == 3) {
					if (rem_pts / (m1 - 1) >= 1) {
						m = m1
					} else {
						m1 = m1 - 1	
					}
				} else {
					m1 = m1 - 1
				}
			}
			#knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
		 	knots = 0:(m-1)/(m-1)*(max(x)-min(x))+min(x)
			#if (pnt) {
			#	knots = 0:(m-1)/(m-1)*(max(xu)-min(xu))+min(xu)
			#} else {
			#	knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#}
		}
	}
	t = 1:(m+4)*0
#new: used for predict
if (is.null(x0)) {
	t[1] = t[2] = min(x)
	t[m+3] = t[m+4] = max(x)
} else {
	t[1] = t[2] = min(x0)
	t[m+3] = t[m+4] = max(x0)
}
	t[3:(m+2)] = knots
##############################
#make quadratic bspline basis#
##############################
	n = length(x)
if (spl) {
	bmat = matrix(0, nrow = n, ncol = (m+1))
#1st edge
	bool = x <= t[4] & x >= t[3]
	bmat[bool, 1] = (t[4] - x[bool]) / (t[4] - t[2]) * (t[4] - x[bool]) / (t[4] - t[3])

#2nd edge 1st part 
	bool = x <= t[4] & x >= t[3]
	bmat[bool, 2] = (x[bool] - t[2]) / (t[4] - t[2]) * (t[4] - x[bool]) / (t[4] - t[3]) + (t[5] - x[bool]) / (t[5] - t[3]) * (x[bool] - t[3]) / (t[4] - t[3]) 

#2nd edge 2nd part 
	bool = x <= t[5] & x >= t[4]
	bmat[bool, 2] = (t[5] - x[bool]) / (t[5] - t[3]) * (t[5] - x[bool]) / (t[5] - t[4])

#3rd edge to the (m-1)th edge
	for(i in 3:(m-1)) {
		#1st part
		bool = x <= t[i+1] & x >= t[i]
		bmat[bool, i] = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	
		#2nd part
		bool = x <= t[i+2] & x >= t[i+1]
		bmat[bool, i] = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])

		#3rd part
		bool = x <= t[i+3] & x >= t[i+2]
		bmat[bool, i] = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	}

#the mth edge 1st part
	bool = x <= t[m+1] & x >= t[m]
	bmat[bool, m] = (x[bool] - t[m])**2 / (t[m+2] - t[m]) / (t[m+1] - t[m])
	 
#the mth edge 2nd part 
	bool = x <= t[m+2] & x >= t[m+1]
	bmat[bool, m] = (x[bool] - t[m]) * (t[m+2] - x[bool]) / (t[m+2] - t[m]) / (t[m+2] - t[m+1]) + (t[m+3] - x[bool]) * (x[bool] - t[m+1]) / (t[m+3] - t[m+1]) / (t[m+2] - t[m+1])

#the (m+1)th edge 
	bool = x <= t[m+2] & x >= t[m+1]
	bmat[bool, m+1] = (x[bool] - t[m+1])**2 / (t[m+3] - t[m+1]) / (t[m+2] - t[m+1])

#####################################################################
#make bqsplines' 1st derivatives at knots. row: knots; column: basis#
#####################################################################
	slopes = matrix(0, nrow = m, ncol = (m+1))
#edge 1
	slopes[1,1] = -2 / (t[4] - t[2])
	#slopes[2,1] = 0 	
	
#edge 2 
	slopes[1,2] = 2 / (t[4] - t[2])
	slopes[2,2] = -2 / (t[5] - t[3])
	#slopes[3,2] = 0

#edge 3 ~ m-1
	for(i in 3:(m-1)){
		#slopes[i-2,i] = 0
		slopes[i-1,i] = 2 / (t[i+2] - t[i])
		slopes[i,i] = -2 / (t[i+3] - t[i+1])
		#slopes[i+1,i] = 0
	}

#edge m
	#slopes[m-2,m] = 0
	slopes[m-1,m] = 2 / (t[m+2] - t[m])
	slopes[m,m] = -2 / (t[m+3] - t[m+1])

#edge m+1
	#slopes[m-1,m+1] = 0
	slopes[m,m+1] = 2 / (t[m+3] - t[m+1])
}
#####################
#plot bqspline basis#
#####################
	xpl = bpl = NULL
if (pic) {
	xpl = 0:1000/1000*(max(x)-min(x))+min(x)
	bpl = matrix(1:(1001*(m+1))*0,nrow=1001)
	
#1st edge
	bool = xpl <= t[4] & xpl >= t[3]
	bpl[bool, 1] = (t[4] - xpl[bool]) / (t[4] - t[2]) * (t[4] - xpl[bool]) / (t[4] - t[3])

#2nd edge 1st part 
	bool = xpl <= t[4] & xpl >= t[3]
	bpl[bool, 2] = (xpl[bool] - t[2]) / (t[4] - t[2]) * (t[4] - xpl[bool]) / (t[4] - t[3]) + (t[5] - xpl[bool]) / (t[5] - t[3]) * (xpl[bool] - t[3]) / (t[4] - t[3]) 

#2nd edge 2nd part 
	bool = xpl <= t[5] & xpl >= t[4]
	bpl[bool, 2] = (t[5] - xpl[bool]) / (t[5] - t[3]) * (t[5] - xpl[bool]) / (t[5] - t[4])

#3rd edge to the (m-1)th edge
	for(i in 3:(m-1)){
		#1st part
		bool = xpl <= t[i+1] & xpl >= t[i]
		bpl[bool, i] = (xpl[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	
		#2nd part
		bool = xpl <= t[i+2] & xpl >= t[i+1]
		bpl[bool, i] = (xpl[bool] - t[i]) * (t[i+2] - xpl[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - xpl[bool]) * (xpl[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])

		#3rd part
		bool = xpl <= t[i+3] & xpl >= t[i+2]
		bpl[bool, i] = (t[i+3] - xpl[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	}

#the mth edge 1st part
	bool = xpl <= t[m+1] & xpl >= t[m]
	bpl[bool, m] = (xpl[bool] - t[m])**2 / (t[m+2] - t[m]) / (t[m+1] - t[m])
	 
#the mth edge 2nd part 

	bool = xpl <= t[m+2] & xpl >= t[m+1]
	bpl[bool, m] = (xpl[bool] - t[m]) * (t[m+2] - xpl[bool]) / (t[m+2] - t[m]) / (t[m+2] - t[m+1]) + (t[m+3] - xpl[bool]) * (xpl[bool] - t[m+1]) / (t[m+3] - t[m+1]) / (t[m+2] - t[m+1])

#the (m+1)th edge 
	bool = xpl <= t[m+2] & xpl >= t[m+1]
	bpl[bool, m+1] = (xpl[bool] - t[m+1])**2 / (t[m+3] - t[m+1]) / (t[m+2] - t[m+1])
}
	ans = new.env()
	ans$bmat = bmat
	ans$bpl = bpl
	ans$xpl = xpl 
	ans$slopes = slopes
	ans$knots = knots
	ans
}


#######################
#bcspline basis matrix#
#######################
bcspl = function(x, m = NULL, knots = NULL, pic = FALSE, pnt = FALSE, spl = TRUE, x0 = NULL) {
############
#make knots#  
############
	bmat = secder = NULL
	n = length(x) #?length(xu)
	xu = unique(x)	
	if (!is.null(knots)) {
		m = length(knots)
	} else {
		if (!is.null(m)) { 
			knots = 0:(m-1)/(m-1)*(max(x)-min(x))+min(x)
			#knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#if (pnt) {
			#	knots = 0:(m-1)/(m-1)*(max(xu)-min(xu))+min(xu)
			#} else {
			#	knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#}
		} else {
			m0 = 4 + round(n^(1 / 9)) 
			m1 = m0
			#m = NULL
			while (is.null(m)) { 
				pts = floor((n - 2) / (m1 - 1))
				rem_pts = (n - 2) %% (m1 - 1)
				if (pts > 3) {
					m = m1
				} else if (pts == 3) {
					if (rem_pts / (m1 - 1) >= 1) {
						m = m1
					} else {
						m1 = m1 - 1	
					}
				} else {
					m1 = m1 - 1
				}
			}
			knots = 0:(m-1)/(m-1)*(max(x)-min(x))+min(x)
			#knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#if (pnt) {
			#	knots = 0:(m-1)/(m-1)*(max(xu)-min(xu))+min(xu)
			#} else {
			#	knots = quantile(xu, probs = seq(0, 1, length = m), names = FALSE)
			#}
		}
	}
	t = 1:(m+6)*0
#new: used for predict
if (is.null(x0)) {
	t[1:3] = min(x)
	t[(m+4):(m+6)] = max(x)
} else {
  t[1:3] = min(x0)
	t[(m+4):(m+6)] = max(x0)
}
	t[4:(m+3)] = knots
##########################
#make cubic bspline basis#
##########################
	#n = length(x)
if (spl) {
	bmat = matrix(0, nrow = n, ncol = m+2)
#1st edge
	i = 1; j = i+1
	bool = x < t[5] & x >= t[4]
	bmat[bool, 1] = (t[5] - x[bool]) / (t[5] - t[2]) * (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	
#2nd edge 1st part
	i = 2; j = i+1
	bool = x <= t[5] & x >= t[4]
	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bmat[bool, 2] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#2nd edge 2nd part
	bool = x <= t[6] & x > t[5]
	b = (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	bmat[bool, 2] = (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 1st part
	i = 3; j = i+1
	bool = x <= t[5] & x >= t[4]
	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bmat[bool, 3] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 2nd part
	bool = x <= t[6] & x > t[5]
	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b =  (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bmat[bool, 3] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 3rd part
	bool = x <= t[7] & x > t[6]
	b = (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	bmat[bool, 3] = (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#4th to (m-1)th edge
	for (i in 4:(m-1)) {
#1st part
		j = i+1
		bool = x <= t[i+1] & x >= t[i]
		a =  (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
		bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a

#2nd part
		bool = x <= t[i+2] & x > t[i+1]
		a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
		b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
		bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#3rd part
		bool = x <= t[i+3] & x > t[i+2]
		a =  (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
		b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
		bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#4th part
		bool = x <= t[i+4] & x > t[i+3]
		b =  (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
		bmat[bool, i] = (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b
}

#mth edge 1st part
	i = m; j = i+1
 	bool = x <= t[m+1] & x >= t[m]
	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a

#mth edge 2nd part
	bool = x <= t[m+2] & x > t[m+1]
	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#mth edge 3rd part
	bool = x <= t[m+3] & x > t[m+2]
	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#(m+1)th edge 1st part
	i = m+1; j = i+1
	bool = x <= t[m+2] & x >= t[m+1]
	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a 

#(m+1)th edge 2nd part
	bool = x <= t[m+3] & x > t[m+2]
	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bmat[bool, i] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b

#(m+2)th edge
	i = m+2
	bool = x <= t[m+3] & x >= t[m+2]
	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bmat[bool, m+2] = (x[bool] - t[i]) / (t[i+3] - t[i]) * a

#####################################
#make bcsplines' 2nd derivatives    #
#at knots. row: knots; column: basis# 
#####################################
	secder = matrix(0, nrow = m, ncol = (m+2))

#1st edge: t4 and t5
	secder[1,1] = 6 / (t[5] - t[2]) / (t[5] - t[3])
	#secder[2,1] = 0 

#2nd edge: t4, t5 and t6
	secder[1,2] = -6 * (1 / (t[5] - t[2]) + 1 / (t[6] - t[3])) / (t[5] - t[3])
	secder[2,2] = 6 / (t[6] - t[3]) / (t[6] - t[4])
	#secder[3,2] = 0 

#3rd edge: t4, t5, t6 and t7
	secder[1,3] = 6 / (t[6] - t[3]) / (t[5] - t[3])	
	secder[2,3] = -6 * (1 / (t[6] - t[3]) + 1 / (t[7] - t[4])) / (t[6] - t[4])
	secder[3,3] = 6 / (t[7] - t[4]) / (t[7] - t[5])
	#secder[4,3] = 0

#4th  to (m-1)th edge: t[i], t[i+1], t[i+2], t[i+3] and t[i+4]
if (m > 4) {
	for (i in 4:(m-1)) {
		#secder[i-3,i] = 0
		secder[i-2,i] = 6 / (t[i+3] - t[i]) / (t[i+2] - t[i])
		secder[i-1,i] = -6 * (1 / (t[i+3] - t[i]) + 1 / (t[i+4] - t[i+1])) / (t[i+3] - t[i+1])
		secder[i,i] = 6 / (t[i+4] - t[i+1]) / (t[i+4] - t[i+2])
		#secder[i+1,i] = 0
	}
}
#mth edge: t[m], t[m+1], t[m+2] and t[m+3]
	#secder[m-3,m] = 0
	secder[m-2,m] = 6 / (t[m+3] - t[m]) / (t[m+2] - t[m])
	secder[m-1,m] = -6 * (1 / (t[m+3] - t[m]) + 1 / (t[m+4] - t[m+1])) / (t[m+3] - t[m+1])
	secder[m,m] = 6 / (t[m+4] - t[m+2]) / (t[m+4] - t[m+1])

#(m+1)th edge: t[m+1], t[m+2] and t[m+3]
	#secder[m-2,m+1] = 0
	secder[m-1,m+1] = 6 / (t[m+4] - t[m+1]) / (t[m+3] - t[m+1])
	secder[m,m+1] = -6 * (1 / (t[m+4] - t[m+1]) + 1 / (t[m+5] - t[m+2])) / (t[m+4] - t[m+2])

#(m+2)th edge: t[m+2] and t[m+3]
	#secder[m-1,m+2] = 0
	secder[m,m+2] = 6 / (t[m+5] - t[m+2]) / (t[m+4] - t[m+2]) 
}
########################
##plot the spline basis# 
########################
xpl = bpl = NULL
if (pic == TRUE) {
	xpl = 0:1000/1000*(max(x)-min(x))+min(x)
	bpl = matrix(0, nrow = 1001, ncol = m+2)
#1st edge
	i = 1; j = i+1
	bool = xpl < t[5] & xpl >= t[4]
	bpl[bool, 1] = (t[5] - xpl[bool]) / (t[5] - t[2]) * (t[j+3] - xpl[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])

#2nd edge 1st part
	i = 2; j = i+1
	bool = xpl <= t[5] & xpl > t[4]
	a = (t[i+3] - xpl[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (xpl[bool] - t[j]) * (t[j+2] - xpl[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - xpl[bool]) * (xpl[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bpl[bool, 2] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#2nd edge 2nd part
	bool = xpl <= t[6] & xpl > t[5]
	b = (t[j+3] - xpl[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	bpl[bool, 2] = (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 1st part
	i = 3; j = i+1
	bool = xpl <= t[5] & xpl > t[4]
	a = (xpl[bool] - t[i]) * (t[i+2] - xpl[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - xpl[bool]) * (xpl[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (xpl[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bpl[bool, 3] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 2nd part
	bool = xpl <= t[6] & xpl > t[5]
	a = (t[i+3] - xpl[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b =  (xpl[bool] - t[j]) * (t[j+2] - xpl[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - xpl[bool]) * (xpl[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bpl[bool, 3] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#3rd edge 3rd part
	bool = xpl <= t[7] & xpl > t[6]
	b = (t[j+3] - xpl[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	bpl[bool, 3] = (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#4th to (m-1)th edge
	for (i in 4:(m-1)) {
		j = i+1
#1st part	
		bool = xpl <= t[i+1] & xpl > t[i]
		a =  (xpl[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
		bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a

#2nd part
		bool = xpl <= t[i+2] & xpl > t[i+1]
		a = (xpl[bool] - t[i]) * (t[i+2] - xpl[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - xpl[bool]) * (xpl[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
		b = (xpl[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
		bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#3rd part
		bool = xpl <= t[i+3] & xpl > t[i+2]
		a =  (t[i+3] - xpl[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
		b = (xpl[bool] - t[j]) * (t[j+2] - xpl[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - xpl[bool]) * (xpl[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
		bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#4th part
		bool = xpl <= t[i+4] & xpl > t[i+3]
		b =  (t[j+3] - xpl[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
		bpl[bool, i] = (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b
}

#mth edge 1st part
	i = m; j = i+1
 	bool = xpl <= t[m+1] & xpl > t[m]
	a = (xpl[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a

#mth edge 2nd part
	bool = xpl <= t[m+2] & xpl > t[m+1]
	a = (xpl[bool] - t[i]) * (t[i+2] - xpl[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - xpl[bool]) * (xpl[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (xpl[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#mth edge 3rd part
	bool = xpl <= t[m+3] & xpl > t[m+2]
	a = (t[i+3] - xpl[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (xpl[bool] - t[j]) * (t[j+2] - xpl[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - xpl[bool]) * (xpl[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#(m+1)th edge 1st part
	i = m+1; j = i+1
	bool = xpl <= t[m+2] & xpl > t[m+1]
	a = (xpl[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a 

#(m+1)th edge 2nd part
	bool = xpl <= t[m+3] & xpl > t[m+2]
	a = (xpl[bool] - t[i]) * (t[i+2] - xpl[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - xpl[bool]) * (xpl[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (xpl[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
	bpl[bool, i] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a + (t[i+4] - xpl[bool]) / (t[i+4] - t[i+1]) * b

#(m+2)th edge
	i = m+2
	bool = xpl <= t[m+3] & xpl > t[m+2]
	a = (xpl[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
	bpl[bool, m+2] = (xpl[bool] - t[i]) / (t[i+3] - t[i]) * a
}
	ans = new.env()
	ans$bmat = bmat
	ans$bpl = bpl
	ans$xpl = xpl
	ans$knots = knots
	ans$secder = secder 
	ans

}

###########################
#bcspline first derivative#
###########################
bcsplfirderi = function(x, knots = NULL, xmin = 0, xmax = 1) {
	m = length(knots) 
	t = 1:(m + 6)*0
	t[1:3] = xmin #min(x)
	t[(m+4):(m+6)] = xmax #max(x)
	t[4:(m+3)] = knots

  	nx = length(x)
	firder =  matrix(0, nrow = nx, ncol = (m + 2)) 

#1st edge: t4 and t5 
	i = 1; j = i+1
	bool = x < t[5] & x >= t[4]
	firder[bool, 1] = -3 * (t[5] - x[bool])^2 / (t[5] - t[2]) / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])

#2nd edge: t4, t5 and t6 
	i = 2; j = i+1
	bool = x <= t[5] & x >= t[4]
	#if (bool) { 
     	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
     	a1 = -2 * (t[i+3] - x[bool]) / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
     	b1 = ((t[j+2] - x[bool])- (x[bool] - t[j])) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + ((t[j+3] - x[bool]) - (x[bool] - t[j+1])) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1]) 
	firder[bool, 2] = a / (t[i+3] - t[i]) + a1 * (x[bool] - t[i]) / (t[i+3] - t[i]) - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
	#}	
	bool = x <= t[6] & x > t[5]
	#if (bool) {
	b = (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
     	b1 =  -2 * (t[j+3] - x[bool]) / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
     	firder[bool, 2] = -b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
	#}

#3rd edge: t4, t5, t6 and t7 
	i = 3; j = i+1
	bool =  x <= t[5] & x >= t[4]
	#if (bool) {
	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j]) 
     	a1 = ((t[i+2] - x[bool]) - (x[bool] - t[i])) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + ((t[i+3] - x[bool]) - (x[bool] - t[i+1])) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
     	b1 = 2 * (x[bool] - t[j]) / (t[j+2] - t[j]) / (t[j+1] - t[j])  
	firder[bool, 3] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1 
	#}	
	bool = x <= t[6] & x > t[5]
	#if (bool) {
	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
     	a1 =  -2 * (t[i+3] - x[bool]) / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
     	b1 = ((t[j+2] - x[bool]) - (x[bool] - t[j])) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + ((t[j+3] - x[bool]) - (x[bool] - t[j+1])) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])   
     	firder[bool, 3] =  a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1 
	#}
	bool = x <= t[7] & x > t[6]
	#if (bool) {
    	b = (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
    	b1 = -2 *(t[j+3] - x[bool]) / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
	firder[bool, 3] = -b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
	#}

#4th  to (m-1)th edge: t[i], t[i+1], t[i+2], t[i+3] and t[i+4] 
if (m > 4) {
	for (i in 4:(m-1)) {
    		j = i+1
		bool = x <= t[i+1] & x >= t[i]
		#if (bool) {
		a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
       		a1 = 2 * (x[bool] - t[i]) / (t[i+2] - t[i]) / (t[i+1] - t[i])
		firder[bool, i] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1
		#}	
		bool = x <= t[i+2] & x > t[i+1]
		#if (bool) {
		a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
		b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
       		a1 =  ((t[i+2] - x[bool]) - (x[bool] - t[i])) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + ((t[i+3] - x[bool]) - (x[bool] - t[i+1])) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
       		b1 = 2 * (x[bool] - t[j]) / (t[j+2] - t[j]) / (t[j+1] - t[j])
		firder[bool, i] =  a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
		#}	
		bool = x <= t[i+3] & x > t[i+2] 
		#if (bool) {
     		a =  (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
 	    	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
      		a1 = -2 * (t[i+3] - x[bool]) / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
      		b1 = ((t[j+2] - x[bool]) - (x[bool] - t[j])) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + ((t[j+3] - x[bool]) - (x[bool] - t[j+1])) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
		firder[bool, i] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
		#}
		bool = x <= t[i+4] & x > t[i+3]
		#if (bool) {
    		b = (t[j+3] - x[bool])**2 / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
      		b1 = -2 * (t[j+3] - x[bool]) / (t[j+3] - t[j+1]) / (t[j+3] - t[j+2])
		firder[bool, i] =  -b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
		#}
	}
}
#mth edge: t[m], t[m+1], t[m+2] and t[m+3] #checked
  	i = m; j = i+1
	bool = x <= t[m+1] & x >= t[m]
	#if (bool) {
    	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
    	a1 = 2 * (x[bool] - t[i]) / (t[i+2] - t[i]) / (t[i+1] - t[i])
	firder[bool, m] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1
	#}	
	bool = x <= t[m+2] & x > t[m+1]
	#if (bool) {
	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
    	a1 =  ((t[i+2] - x[bool]) - (x[bool] - t[i])) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + ((t[i+3] - x[bool]) - (x[bool] - t[i+1]) ) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1]) 
    	b1 = 2 * (x[bool] - t[j]) / (t[j+2] - t[j]) / (t[j+1] - t[j])
	firder[bool, m] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1
	#}	
	bool = x <= t[m+3] & x > t[m+2] 
	#if (bool) {
    	a = (t[i+3] - x[bool])**2 / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
	b = (x[bool] - t[j]) * (t[j+2] - x[bool]) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + (t[j+3] - x[bool]) * (x[bool] - t[j+1]) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
    	a1 =  -2 * (t[i+3] - x[bool]) / (t[i+3] - t[i+1]) / (t[i+3] - t[i+2])
    	b1 = ((t[j+2] - x[bool]) - (x[bool] - t[j])) / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) + ((t[j+3] - x[bool]) - (x[bool] - t[j+1])) / (t[j+3] - t[j+1]) / (t[j+2] - t[j+1])
	firder[bool, m] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1  - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1  
	#}	

#(m+1)th edge: t[m+1], t[m+2] and t[m+3] 
  	i = m+1; j = i+1
	bool = x <= t[m+2] & x >= t[m+1]
	#if (bool) {
    	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
    	a1 = 2 * (x[bool] - t[i]) / (t[i+2] - t[i]) / (t[i+1] - t[i])
	firder[bool, m+1] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1
	#}
	bool = x <= t[m+3] & x > t[m+2]
	#if (bool) {
    	a = (x[bool] - t[i]) * (t[i+2] - x[bool]) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + (t[i+3] - x[bool]) * (x[bool] - t[i+1]) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
	b = (x[bool] - t[j])**2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
    	a1 = ((t[i+2] - x[bool]) - (x[bool] - t[i])) / (t[i+2] - t[i]) / (t[i+2] - t[i+1]) + ((t[i+3] - x[bool]) - (x[bool] - t[i+1])) / (t[i+3] - t[i+1]) / (t[i+2] - t[i+1])
    	b1 = 2 * (x[bool] - t[j]) / (t[j+2] - t[j]) / (t[j+1] - t[j])
	firder[bool, m+1] = a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1 - b / (t[i+4] - t[i+1]) + (t[i+4] - x[bool]) / (t[i+4] - t[i+1]) * b1  
	#}

#(m+2)th edge: t[m+2] and t[m+3]
  	i = m+2
	bool = x <= t[m+3] & x >= t[m+2]
	#if (bool) {
    	a = (x[bool] - t[i])**2 / (t[i+2] - t[i]) / (t[i+1] - t[i])
    	a1 = 2 * (x[bool] - t[i]) / (t[i+2] - t[i]) / (t[i+1] - t[i]) 
	firder[bool, m+2] =  a / (t[i+3] - t[i]) + (x[bool] - t[i]) / (t[i+3] - t[i]) * a1
	#}	
	firder
}

###########
#CicFamily#
###########
CicFamily <- function(object,...)UseMethod("CicFamily")
CicFamily <- function(object) {
  llh.fun <- function(y, muhat = NULL, etahat = 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
      }
    }
    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 
    }
    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])}
           }         
         }
    }
    gr
  }

  wt.fun <- function(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 
    }
    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])
    }
    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
      #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
    }
   muhat 
  }

  ysim.fun <- function(n, mu0 = NULL, fml = object$family) {
    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)
    }
    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
    }
   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)
  }
###################
#get null deviance#
###################
  if (fml == "binomial" | fml == "poisson") {
      diff <- 1
      muhat0 <- mean(y) + 1:n*0
      if (fml == "poisson") {
         etahat0 <- log(muhat0)
      } 
      if (fml == "binomial") {
         etahat0 <- log(muhat0 / (1 - 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 == "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)
  class(ans) <- "CicFamily"
  return (ans)
}

##########
#make_arp#
##########
make_arp = function(y, pnt = FALSE, dmat = NULL, bmat = NULL, sl = NULL, p_max = 2, m2 = NULL, lambdas = NULL, pen_fit = NULL, p_fit = NULL, dv0 = NULL, qv0 = NULL, cv = NULL, hs = NULL, ids = NULL, constr = TRUE) {
	if (is.null(pen_fit) & is.null(p_fit)) {
		ps = 1:p_max
		if (!pnt) {
	   		lambdas = as.vector(0)
 		}
  		la = length(lambdas)
#include p = 0
  		nposs = la * (p_max + 1)
	} else {
  		ps = as.vector(p_fit)
  		lambdas = as.vector(pen_fit)
  		la = p_max = nposs = 1
	}
  	aics = aics1 = aics2 = trs = trrs = sigs = pens = 1:nposs*0
	n = length(y)
  	phis = bhs = r_mats = pmat0s = list()
	fhats = fhats_x = ess = rss = matrix(0, nrow = nposs, ncol = n)
	sm = 1e-5
	iter = 0
# for p: 0 ~ 4
  	rs = rs2 = 1:n*0
  	r_mat = tau_mat = matrix(0, n, n)
	qv0 = crossprod(bmat)
	tdf = ncol(bmat)
    	imat = diag(tdf)
  	for (ila in 1:la) {
        	pen = lambdas[ila]
         	if (pnt) {
            		qv = qv0 + pen * dv0
         	} else {
            		qv = qv0
         	}	
		cv = crossprod(bmat, y)
#new, for jp:
		if (constr) {
         		qans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)	
			bh = fitted(qans)
         		qdf = qans$df
          		fhat0 = bmat %*% bh
          		d_use = t(qans$xmat)
		} else {
			bh = solve(qv0, cv)
			qdf = tdf 
			fhat0 = bmat %*% bh	
			d_use = bmat
		}          	
		if (qdf < tdf) {
             		pd = d_use %*% solve(crossprod(d_use), t(d_use))
             		pmat0 = imat - pd
          	} else {
            		pmat0 = imat
          	}                    
          	umat = chol(qv)
          	iumat = diag(ncol(umat))
          	uinv = backsolve(umat, iumat)
          	bu = bmat %*% uinv    
          	pmat = bu %*% tcrossprod(pmat0, bu)
#constrained
          	tr = sum(diag(pmat))
		pu = bmat %*% solve(qv, t(bmat))
#unconstrained
		trr = sum(diag(pu))  
          	es = y - fhat0
          	nd = ncol(d_use)
          	sig2 =  sum((y - fhat0)^2) / (n - 1.5 * (tdf - nd))
          	aic = n * log(sig2) + 2 * tr
  		if (is.null(p_fit) || p_fit == 0) {
   			iter = iter + 1
        		phis[[iter]] = 0
         		bhs[[iter]] = bh
			r_mats[[iter]] = diag(n)    		
			fhats[iter, ] = fhat0
         		ess[iter, ] = es
         		aics[iter] = aic		
         		trs[iter] = tr
			trrs[iter] = trr 
         		sigs[iter] = sig2
         		pens[iter] = pen    
		}
		if (is.null(p_fit) || p_fit != 0) {
			for (ip in 1:p_max) {
				nrep = 0
  				sm0 = 1
                		fhat = fhat0 
				iter = iter + 1  				
				p = ps[ip]  
  				while (sm0 > sm & nrep < 1e+3) {
                      			nrep = nrep + 1
                      			oldf = fhat 
                      			es = y - fhat
					es = es - mean(es)                      			
					rs = sapply(hs, FUN = function(h, X = es, N = n) {sum((X[1:(N - h)] - mean(X)) * (X[(1 + h):N]- mean(X))) / N})
					tau_mat = matrix(0, n, n)
                      			for (ih in 1:n) {
                        		  	ids_i = ids[[ih]]
	                        		h = hs[ih]
	                        		tau_mat[ids_i] = rs[h + 1]
                      			}               
					tau_p_mat = tau_mat[1:p, 1:p]
            				rp = rs[2:(p + 1)]
					phi = solve(tau_p_mat, rp)  
					sig2 = rs[1] - crossprod(phi, rp)        
					rs2 = make_rs(phi, ord = p, nvec = n, sig2 = sig2)
					r_mat = matrix(0, n, n)
     					for (ih in 1:n) {
                        		  	ids_i = ids[[ih]]
	                        		h = hs[ih]
	                        		r_mat[ids_i] = rs2[h + 1]
                      			}
#this line makes the fit flat!
					r_mat = r_mat / rs2[1]
					if (pnt) {
     	                   			qv = crossprod(bmat, solve(r_mat, bmat)) + pen * dv0
    		        		} else {
         					qv = crossprod(bmat, solve(r_mat, bmat))                        	
      		        		}
					cv = crossprod(bmat, solve(r_mat, y))
#new, for jp:
					if (constr) {
                      				qans = qprog(qv, cv, sl, 1:nrow(sl)*0, msg = FALSE)
                      				bh = fitted(qans)
                      				fhat = bmat %*% bh
					} else {
						bh = solve(qv, cv)
						fhat = bmat %*% bh	
					}
#lines(x, fhat, lty = 2, col = nrep)
                      			sm0 = mean((fhat - oldf)^2)	
            			}
				if (constr) {		
					d_use = t(qans$xmat)
            				qdf = qans$df
				} else {
					d_use = bmat
					qdf = tdf 
				}
            			if (qdf < tdf) {
					pd = d_use %*% solve(crossprod(d_use), t(d_use))
               				pmat0 = imat - pd
            			} else {
              				pmat0 = imat
            			}   
				umat = chol(qv)
				iumat = diag(ncol(umat))
            			uinv = backsolve(umat, iumat)
           			bu = bmat %*% uinv
            			pmat = bu %*% tcrossprod(pmat0, bu) %*% chol2inv(chol(r_mat))    
				pu = bmat %*% solve(qv, crossprod(bmat, chol2inv(chol(r_mat))))
#constrained
            			tr = sum(diag(pmat))
#unconstrained	
				trr = sum(diag(pu))   			
				rss[iter, ] = rs2
				pmat0s[[iter]] = pmat0         		
         			bhs[[iter]] = bh
	 			fhats[iter, ] = fhat              
         			ess[iter, ] = es 
        			r_mats[[iter]] = r_mat
        			phis[[iter]] = phi
	 			aics1[iter] = n * log(sig2)
            			aics2[iter] = 2 * (p + tr) 	
            			aics[iter] = n * log(sig2) + 2 * (p + tr)	
         			trs[iter] = tr
				trrs[iter] = trr
         			sigs[iter] = sig2
         			pens[iter] = pen
   			}
     		}                     
	}
  	id_aic = which(aics == min(aics))
  	lid = length(id_aic)
  	if (lid > 1) {
     		id_use = id_aic[1]
  	} else {
     		id_use = id_aic
  	}
#print (id_use)
	phi_use = phis[[id_use]]
	r_use = r_mats[[id_use]]
  	bh_use = bhs[[id_use]]
	fhat_use = fhats[id_use, ]
  	tr_use = trs[id_use]
	trr_use = trrs[id_use]
  	sig_use = sigs[id_use]
  	es_use = ess[id_use, ]
  	lambda_use = pens[id_use]
#lambda by order
	sigs = matrix(sigs, nrow = la, byrow = TRUE)
	aics = matrix(aics, nrow = la, byrow = TRUE)
	trs = matrix(trs, nrow = la, byrow = TRUE)
	trrs = matrix(trrs, nrow = la, byrow = TRUE)
	pens = matrix(pens, nrow = la, byrow = TRUE) 	
	aicmin = min(aics)
  	rslt = list(phi_use = phi_use, phis = phis, thetahat = bh_use, fhat_use = fhat_use, fhats = fhats, df = tr_use, r_use = r_use, sig_use = sig_use, sigs = sigs, pens = pens, trs = trs, trrs = trrs, aics = aics, aics1 = aics1, aics2 = aics2, aicmin = aicmin, tr_use = tr_use, trr_use = trr_use, es_use = es_use, ess = ess, lambda_use = lambda_use)
  	rslt
}


############################
#make theoretical residuals#
############################
make_rs = function(phi, ord = 1, nvec, sig2) {
	rhs = rvec = 1:nvec*0
	if (ord == 1) {
		r0 = sig2 / (1 - phi^2)
		rhs[1] = 1
		rhs[2:nvec] = phi^(1:(nvec - 1))
		rvec = rhs * r0	
		#rvec = rhs	
	}
	if (ord == 2) {
		phi1 = phi[1]; phi2 = phi[2]
		rh1 = phi1 / (1 - phi2); rh2 = phi1^2 / (1 - phi2) + phi2
		rhs[1] = 1; rhs[2] = rh1; rhs[3] = rh2
		rh_vec = c(rh1, rh2)
		r0 = sig2 / (1 - sum(phi * rh_vec))
		if (nvec >= 3) {
			rh_iter = rev(rh_vec)
			for (l in 1:(nvec - 3)) {
				rhs[3 + l] = sum(phi * rh_iter)
				rh_iter = c(rhs[3 + l], rh_iter[1])
				if (all(abs(rh_iter) < 1e-10)) {			
					break
				}
			}	
		}
		rvec = rhs * r0
		#rvec = rhs	
	}
	if (ord == 3) {
		phi_mat = matrix(0, 3, 3)
		phi_mat[1, 1] = 1 - phi[2]; phi_mat[1, 2] = -phi[3]
		phi_mat[2, 1] = -phi[1] - phi[3]; phi_mat[2, 2] = 1
		phi_mat[3, 1] = -phi[2]; phi_mat[3, 2] = -phi[1]; phi_mat[3, 3] = 1
		rh_vec = solve(phi_mat, phi)
		rhs[1] = 1
		rhs[2:4] = rh_vec
		r0 = sig2 / (1 - sum(phi * rh_vec))
		if (nvec >= 4) {
			rh_iter = rev(rh_vec)
			for (l in 1:(nvec - 4)) {
				rhs[4 + l] = sum(phi * rh_iter)
				rh_iter = c(rhs[4 + l], rh_iter[1:2])
				if (all(abs(rh_iter) < 1e-10)) {
					break
				}
			}	
		}
		rvec = rhs * r0
		#rvec = rhs	
	}
	if (ord == 4) {
		phi_mat = matrix(0, 4, 4)
		phi_mat[1, 1] = 1 - phi[2]; phi_mat[1, 2] = -phi[3]; phi_mat[1, 3] = -phi[4]
		phi_mat[2, 1] = -phi[1] - phi[3]; phi_mat[2, 2] = -phi[4] + 1
		phi_mat[3, 1] = -phi[2] - phi[4]; phi_mat[3, 2] = -phi[1]; phi_mat[3, 3] = 1
		phi_mat[4, 1] = -phi[3]; phi_mat[4, 2] = -phi[2]; phi_mat[4, 3] = -phi[1]; phi_mat[4, 4] = 1
		rh_vec = solve(phi_mat, phi)
		rhs[1] = 1
		rhs[2:5] = rh_vec
		r0 = sig2 / (1 - sum(phi * rh_vec))
		if (nvec >= 5) {
			rh_iter = rev(rh_vec)
			for (l in 1:(nvec - 5)) {
				rhs[5 + l] = sum(phi * rh_iter)
				rh_iter = c(rhs[5 + l], rh_iter[1:3]) 
				if (all(abs(rh_iter) < 1e-10)) {
					break
				}
			}	
		}
		rvec = rhs * r0
		#rvec = rhs
	}
	rvec = round(rvec, 8)
	return (rvec)
}

###
make_dmat = function(m, q, pnt = FALSE) {
	dmat = NULL
	if (pnt) {
		dmat = matrix(0, nrow = (m - q), ncol = m)
#  third-order
		if (q == 3) {
			for (i in 4:m) {
				dmat[i-3, i-3] = 1; dmat[i-3, i-2] = -3; dmat[i-3, i-1] = 3; dmat[i-3, i] = -1
			}
		}
# second order
		if (q == 2) {
			for (i in 3:m) {
				dmat[i-2, i-2] = 1; dmat[i-2, i-1] = -2; dmat[i-2, i] = 1
			}
		}
# first order
		if (q == 1) {
			for (i in 2:m) {
				dmat[i-1, i-1] = 1; dmat[i-1, i] = -1
			}
		}
# zero order
		if (q == 0) {
			for (i in 1:m) {
				dmat[i, i] = 1
			}
		}
	}
	return (dmat)
}


####
find_chpt = function(dist, pos, m2, knots) {
  #print (dist)
  sub = NULL
  sm = 1e-5
  if (abs(dist[pos]) < sm) {
    if (pos == 1) {
      obs.r = 2:m2
      dist.r = dist[obs.r]
      if (dist.r[1] < 0) {
        x.l = knots[pos]
        x.r = knots[pos]
      } else if (abs(dist.r[1]) < sm) {
        if (any(dist.r < 0)) {
          id.r = (obs.r[dist.r < 0])[1] - 1
          x.l = knots[pos]
          x.r = knots[id.r]
        } else {
          x.l = knots[pos]
          x.r = knots[m2]
        }
      }
    } else if (pos == m2) {
      obs.l = 1:(m2 - 1)
      dist.l = dist[obs.l]
      if (dist.l[m2 - 1] > sm) {
        x.l = knots[pos]
        x.r = knots[pos]
      } else if (abs(dist.l[m2 - 1]) < sm) {
        if (any(dist.l > sm)) {
          id.l = tail(obs.l[dist.l > sm], 1) + 1
          x.l = knots[id.l] 
          x.r = knots[pos]
        } else {
          x.l = knots[1]
          x.r = knots[pos]
        }
      }
    } else {
      obs.l = 1:(pos - 1) 
      obs.r = (pos + 1):m2
      dist.l = dist[obs.l]
      dist.r = dist[obs.r]
      if (dist.l[pos - 1] > sm) {
        x.l = knots[pos]
      } else if (abs(dist.l[pos - 1]) < sm) {
        if (any(dist.l > sm)) {
          id.l = tail(obs.l[dist.l > sm], 1) + 1
          x.l = knots[id.l] 
          x.r = knots[pos]
        } else {
          x.l = knots[1] 
          x.r = knots[pos]
        }	
      }
      if (dist.r[1] < 0) {
        x.r = knots[pos]
      } else if (abs(dist.r[1]) < sm) {
        if (any(dist.r < (-sm))) {
          id.r = (obs.r[dist.r < 0])[1] - 1
          x.l = knots[pos] 
          x.r = knots[id.r]
        } else {
          x.l = knots[pos]
          x.r = knots[m2]
        }
      }
    }
    chpt = (x.l + x.r) / 2
  } 
  if (dist[pos] < (-sm)) {
    if (pos == 1) {
      chpt = knots[pos]
    } else {
      if (any(dist > sm)) {
        if (dist[pos - 1] > sm) {
          x.l = knots[pos - 1] 
          x.r = knots[pos]
          y.l = dist[pos - 1] 
          y.r = dist[pos]
          sub = c(x.l, x.r)
          a = (y.r - y.l) / (x.r - x.l)
          b = y.r - a * x.r 
          #chpt = - (y.r - ((y.r - y.l) / (x.r - x.l)) * x.r ) / (y.r - y.l) / (x.r - x.l)
          f = function(x) a * x + b
          chpt = uniroot(f, c(x.l, x.r), tol = 1e-8)$root
        } 
        if (abs(dist[pos - 1]) < sm) {
          obs.l = 1:(pos - 1)
          dist.l = dist[obs.l]
          id.l = obs.l[dist.l > sm]
          x.l = knots[tail(id.l, 1) + 1]
          x.r = knots[pos - 1]
          chpt = (x.l + x.r) / 2
          #sub = c(x.l, x.r)
          y.l = dist[tail(id.l, 1) + 1]
          y.r = dist[pos - 1]
        }
      } else if (all(abs(dist[1:(pos - 1)]) < sm)) {
        x.l = knots[1]
        x.r = knots[pos - 1]
        chpt = (x.l + x.r) / 2
        y.l = dist[1]
        y.r = dist[pos]
        #sub = c(x.l, x.r)
      }
    }
  }
  return (chpt)	
}

Try the ShapeChange package in your browser

Any scripts or data that you put into this service are public.

ShapeChange documentation built on Sept. 8, 2023, 6:06 p.m.