.packageName <- "sampling"
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/HTestimator.R"
HTestimator<-function(y,pik) 
{
if(any(is.na(pik))) stop("there are missing values in pik")
if(any(is.na(y))) stop("there are missing values in y")
if(length(y)!=length(pik)) stop("y and pik have different sizes")
crossprod(y,1/pik)
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/HTstrata.r"
HTstrata<-function (y, pik, strata, description=FALSE) 
{
    str <- function(st, h, n) .C("str", as.double(st), as.integer(h), 
        as.integer(n), s = double(n), PACKAGE = "sampling")$s
    if (is.matrix(y)) 
        sample.size <- nrow(y)
    else sample.size <- length(y)
    h <- unique(strata)
    s1 <- 0
    for (i in 1:length(h)) {
        s <- str(strata, h[i], sample.size)
        est<-HTestimator(y[s == 1], pik[s == 1])
        s1 <- s1 + est
        if(description)
         cat("For stratum ",i,", the Horvitz-Thompson estimator is:",est,"\n")
    }
if(description)
  cat("The Horvitz-Thompson estimator is:\n")
s1
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPMEpik2frompikw.R"
"UPMEpik2frompikw" <-
function(pik,w)
{
n=sum(pik)
N=length(pik)
M=array(0,c(N,N))
for(k in 1:N) for(l in 1:N)
  if(pik[k]!=pik[l] & k!=l) M[k,l]= (pik[k]*w[l]-pik[l]*w[k])/(w[l]-w[k]) else M[k,l]=-1
for(i in 1:N) M[i,i]=pik[i]
for(k in 1:N)
  {
  tt=0
  comp=0
  for(l in 1:N)
        {if(M[k,l]!=-1) tt=tt+M[k,l]
         else comp=comp+1;
        }
        cc=(n*pik[k]-tt)/comp
for(l in 1:N)  if(M[k,l]==-1) M[k,l]=cc
  }
M
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPMEpikfromq.R"
"UPMEpikfromq" <-
function(q) 
{ 
n=ncol(q)
N=nrow(q)
pro=array(0,c(N,n))
pro[1,n]=1 
for(i in 2:N) 
for(j in 2:n) 
{ 
pro[i,j]=pro[i,j]+pro[i-1,j]*(1-q[i-1,j]) 
pro[i,j-1]=pro[i,j-1]+pro[i-1,j]*(q[i-1,j]) 
}; 
for(i in 2:N) 
{ 
pro[i,1]=pro[i,1]+pro[i-1,1]*(1-q[i-1,1]) 
} 
rowSums(pro*q) 
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPMEpiktildefrompik.R"
"UPMEpiktildefrompik" <-
function(pik,eps=1e-6) 
{ 
n=sum(pik) 
pikt=pik 
arr=1 
while(arr>eps) 
{ 
w=(pikt)/(1-pikt) 
q=UPMEqfromw(w,n)
pikt1=pikt+pik-UPMEpikfromq(q) 
arr=sum(abs(pikt-pikt1)) 
pikt=pikt1 
} 
pikt
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPMEqfromw.R"
"UPMEqfromw" <-
function(w,n)
{
N=length(w)
expa=array(0,c(N,n))
for(i in 1:N) expa[i,1]= sum(w[i:N]) 
for(i in (N-n+1):N) expa[i,N-i+1]=exp(sum(log(w[i:N]))) 
for(i in (N-2):1) 
for(z in 2:min(N-i,n)) 
{ 
expa[i,z]=w[i]*expa[i+1,z-1]+expa[i+1,z] 
} 
q=array(0,c(N,n))
for(i in N:1) q[i,1]= w[i]/expa[i,1]
for(i in N:(N-n+1)) q[i,N-i+1]=1
for(i in (N-2):1) 
for(z in 2:min(N-i,n)) 
q[i,z] = w[i]*expa[i+1,z-1]/expa[i,z]
q
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPMEsfromq.R"
"UPMEsfromq" <-
function(q)
{
n=ncol(q)
N=nrow(q)
s=rep(0,times=N) 
for(k in 1:N) 
if(n!=0) if(runif(1)<q[k,n]) {s[k]=1;n=n-1} 
s 
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPbrewer.R"
"UPbrewer" <- function(pik, eps = 1e-06)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
list = pik > eps & pik < 1 - eps
pikb = pik[list]
N = length(pikb)
s=pik
if(N<1) stop("the pik vector has all elements outside of the range [eps,1-eps]")
else 
{
sb=rep(0,N)
n=sum(pikb)
for (i in 1:n) {
        a = sum(pikb*sb)
        p = (1-sb)*pikb*((n-a)-pikb)/((n-a)-pikb*(n-i+1))
        p = p/sum(p)
        p = cumsum(p)
        u=runif(1)
        for(j in 1:length(p))
		if(u<p[j]) break;
        sb[j] = 1
}
s[list]=sb
}
s
}



#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPmaxentropy.R"
"UPmaxentropy" <-
function(pik) 
{ 
n=sum(pik)
n=.as_int(n)
if(n>=2)
{
pik2=pik[pik!=1]
n=sum(pik2)
piktilde=UPMEpiktildefrompik(pik2) 
w=piktilde/(1-piktilde) 
q=UPMEqfromw(w,n)
s2=UPMEsfromq(q)
s=rep(0,times=length(pik))
s[pik==1]=1
s[pik!=1][s2==1]=1
}
if(n==0) s=rep(0,times=length(pik))
if(n==1) s=as.vector(rmultinom(1, 1,pik)) 
s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPmaxentropypi2.R"
"UPmaxentropypi2" <-
function(pik)
{
n=sum(pik)
N=length(pik)
M=array(0,c(N,N))
if(n>=2)
{
pik2=pik[pik>0 & pik<1]
pikt=UPMEpiktildefrompik(pik2)
w=pikt/(1-pikt)
M[pik>0 & pik<1,pik>0 & pik<1]=UPMEpik2frompikw(pik2,w)
M[,pik==1]=pik
for(k in 1:N) if(pik[k]==1)  M[k,]=pik
}
if(n==1) for(k in 1:N) M[k,k]=pik[k];
M
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPmidzuno.R"
"UPmidzuno" <-
function(pik,eps=1e-6) 1-UPtille(1-pik,eps)

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPmidzunopi2.R"
"UPmidzunopi2" <-
function(pik) 
{
N=length(pik)
UN=rep(1,times=N)
b=1-pik%*%t(UN)
1-b-t(b)+UPtillepi2(1-pik)
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPminimalsupport.R"
"UPminimalsupport" <-
function(pik)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
basicsplit<-function(pik)
{
N=length(pik)
n=sum(pik)
A=(1:N)[pik==0]
B=(1:N)[pik==1]
C=setdiff(setdiff(1:N,A),B)
D=C[sample(length(C), round(n-length(B)))]
s1v=rep(0,times=N)
s1v[c(B,D)]=1
alpha=min(1-max(pik[setdiff(C,D)]),min(pik[D]))
pikb= (pik-alpha*s1v)/(1-alpha)
if(runif(1,0,1)<alpha) s=s1v else s=pikb
s
}
is.a.sample<-function(s,EPS=sqrt(.Machine$double.eps)) if(sum(abs(s-round(s)))<EPS) TRUE else FALSE
while(!is.a.sample(pik))pik=basicsplit(pik) 
round(pik)
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPmultinomial.R"
"UPmultinomial" <-
function(pik) 
{if(any(is.na(pik))) stop("there are missing values in the pik vector")
as.vector(rmultinom(1,sum(pik),pik/sum(pik)))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPopips.r"
UPopips<-function(lambda, type=c("pareto","uniform","exponential"))
{
if(any(is.na(lambda))) stop("there are missing values in the lambda vector")
n=sum(lambda)
if(!(type %in% c("pareto","uniform","exponential"))) stop("the type argument is wrong")
omega=runif(n)
switch(type,pareto=order(omega*(1-lambda)/((1-omega)*lambda))[1:n],
uniform=order(omega/lambda)[1:n], exponential=order(log(1-omega)/log(1-lambda)))[1:n]
}




#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPpivotal.r"
"UPpivotal" <-
function(pik,eps=1e-6)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
N<-length(pik)
s<-rep(0,times=N)
a<-pik[1]
b<-pik[2]
i<-1
j<-2
k<-3
while(k<=N)
{
u<-runif(1)
if(a>=eps & a<= 1-eps & b>=eps & b<= 1-eps)
if(a+b>1)
	{
		if(u<(1-b)/(2-a-b)) 
			{b<-a+b-1;a<-1} 
		else {a<-a+b-1;b<-1} 
      }
	else{ if(u< b/(a+b)) 
			{b<- a+b;a<-0} 
            else {a<- a+b;b<-0} 
           }
      if( (a<eps | a > 1-eps)& (k<=N)) 
                {s[i]=a;a=pik[k];i=k;k=k+1;} 
      if( (b<eps | b > 1-eps)& (k<=N) ) 
                  {s[j]=b;b=pik[k];j=k;k=k+1;} 
}
u<-runif(1)
if(a>=eps & a<= 1-eps & b>=eps & b<= 1-eps)
if(a+b>1)
          {
            if(u<(1-b)/(2-a-b)) {b<-a+b-1;a<-1} 
            else {a<-a+b-1;b<-1} 
          }
 else{ if(u< b/(a+b)) 
                 		{b<-a+b;a<-0} 
       else {a<- a+b;b<-0} 
       } 
s[i]=a; s[j]=b;
s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPpoisson.R"
"UPpoisson" <-
function(pik) 
{if(any(is.na(pik))) stop("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
as.numeric(runif(length(pik))<pik)
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPrandompivotal.R"
"UPrandompivotal" <-
function(pik,eps=1e-6)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
N=length(pik)
v=sample(c(N),N)
s=numeric(N)
s[v]=UPpivotal(pik[v],eps)
s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPrandomsystematic.R"
"UPrandomsystematic" <-
function(pik,eps=1e-6)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
N=length(pik)
v=sample(N,N)
s=numeric(N)
s[v]=UPsystematic(pik[v],eps)
s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPsampford.R"
"UPsampford" <-
function(pik,eps=1e-6)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
list= pik>eps & pik < 1-eps
pikb=pik[list]
n=sum(pikb)
N=length(pikb)
s=pik
if(N<1) stop("the pik vector has all elements outside of the range [eps,1-eps]")
else 
{
sb=rep(2,N)
y=pikb/(1-pikb)/sum(pikb/(1-pikb))
while(sum(sb<=1)!=N)
	sb=as.vector(rmultinom(1,1,pikb/sum(pikb))+rmultinom(1,n-1,y))
s[list]=sb
}
s
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPsampfordpi2.r"
UPsampfordpi2=function(pik) 
{
n=sum(pik)
n=.as_int(n)
N=length(pik)
p=pik/n
pikl=matrix(0,N,N)
Lm=rep(0, n)
lambda=p/(1-n*p)
Lm[1]=1
if(n>=2)
for (i in 2:n) {
for (r in 1:(i-1)) Lm[i]=Lm[i]+((-1)^(r-1))*sum(lambda^r)*Lm[i-r]
Lm[i]=Lm[i]/(i - 1)
}
t1=(n + 1) - (1:n)
Kn=1/sum(t1*Lm/n^t1)
Lm2=rep(0, n - 1)
t2=(1:(n - 1))
t3=n - t2 
for (i in 2:N) {
for (j in 1:(i - 1)) {
Lm2[1]=1
Lm2[2]=Lm[2] - (lambda[i] + lambda[j])
if(n>3)
for (m in 3:(n - 1)) {
Lm2[m]=Lm[m] - (lambda[i] + lambda[j]) * Lm2[m -1] - lambda[i] * lambda[j] * Lm2[m - 2]
}
pikl[i, j]=Kn * lambda[i] * lambda[j] * sum((t2+1-n*(p[i] + p[j]))*Lm2[t3]/n^(t2 - 1))
pikl[j, i]=pikl[i, j]
}
pikl[i, i]=pik[i] 
}
pikl[1, 1]=pik[1] 
pikl
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPsystematic.R"
"UPsystematic"<-function(pik,eps=1e-6) 
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
list=pik > eps & pik < 1-eps
pik1 = pik[list]
N = length(pik1)
a = (c(0, cumsum(pik1)) - runif(1, 0, 1))%%1
s1 = as.integer(a[1:N] > a[2:(N + 1)])
s = pik
s[list] = s1
s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPsystematicpi2.R"
UPsystematicpi2<-function(pik) 
{
    n=sum(pik)
    n=.as_int(n)
    pik1 = pik[pik > 0 & pik < 1]
    N = length(pik1)
    Vk = cumsum(pik1)
    Vk1=Vk%%1 
    if(Vk1[N]!=0) Vk1[N]=0
    r = c(sort(Vk1), 1)
    cent = (r[1:N] + r[2:(N + 1)])/2
    p = r[2:(N + 1)] - r[1:N]
    A = matrix(c(0, Vk), nrow = N + 1, ncol = N) - 
    t(matrix(cent,nrow = N, ncol = N + 1))
    A = A%%1
    M = matrix(as.integer(A[1:N, ] > A[2:(N + 1), ]), N, N)
    pi21 = M %*% diag(p) %*% t(M)
    pi2 = pik %*% t(pik)
    pi2[pik > 0 & pik < 1, pik > 0 & pik < 1] = pi21
    pi2
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPtille.R"
"UPtille" <-
function(pik,eps=1e-6)
{
if(any(is.na(pik))) stop("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
list = pik > eps & pik < 1 - eps
pikb = pik[list]
N = length(pikb)
s=pik
if(N<1) stop("the pik vector has all elements outside of the range [eps,1-eps]")
else 
{
n=sum(pikb)
sb=rep(1,N)
b=rep(1,N)
for(i in 1:(N-n))
	{a=inclusionprobabilities(pikb,N-i)
	v=1-a/b
	b=a
	p=v*sb
	p=cumsum(p)
	u=runif(1)
        for(j in 1:length(p))
		if(u<p[j]) break
        sb[j] = 0
	}
s[list]=sb
}
s
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/UPtillepi2.R"
"UPtillepi2" <-
function(pik,eps=1e-6)
{

if(any(is.na(pik))) warning("there are missing values in the pik vector")
n=sum(pik)
n=.as_int(n)
list = pik > eps & pik < 1 - eps
pikb = pik[list]
N = length(pikb)
#ppf=pik%*%t(pik)
ppf=matrix(0,length(pik),length(pik))
if(N<1) stop("the pik vector has all elements outside of the range [eps,1-eps]")
else 
{
n=sum(pikb)
if(N>n)
{
UN=rep(1,N)
b=rep(1,N)
pp=1
for(i in 1:(N-n))
	{
	a=inclusionprobabilities(pikb,N-i)
	vv=1-a/b
	b=a
	d=vv %*% t(UN)
	pp=pp*(1-d-t(d))
	}
diag(pp)=pikb
ppf[list,list]=pp
}
}
ppf
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/as_int.r"
.as_int<-function(x)
{if(!is.integer(x)) {
            xo =x
            x = round(x)
            if(any(x > .Machine$integer.max))
                stop("the input has entries too large to be integer")
            if(!identical(TRUE, (ax <- all.equal(xo, x))))
                warning("the sum of pik has been rounded to integer: ", ax)
            storage.mode(x) <- "integer"
        }
x
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/balancedcluster.R"
"balancedcluster" <-
function(X,m,cluster,selection=1,comment=TRUE,method=1)
{
cluster=cleanstrata(cluster)
if(comment==TRUE) cat("\nSELECTION OF A SAMPLE OF CLUSTERS\n")
p=dim(X)[2]
N=dim(X)[1]
H=max(cluster)
XC=array(0,c(H,p))
Ni=rep(0,times=H)
for(h in 1:H)
   { Ni[h]=sum(as.integer(cluster==h)) 
    for(j in 1:p)  XC[h,j]=sum(X[cluster==h,j]) 
   }
if(selection==1) pik=inclusionprobabilities(Ni,m) else pik=rep(m/H,times=H)
s=samplecube(cbind(pik,XC),pik,1,comment,method) 
res=array(0,c(N,2))
for(h in 1:H) 
    {
     res[cluster==h,1]=s[h] 
     res[cluster==h,2]=pik[h]
     }
res
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/balancedstratification.R"
"balancedstratification" <-
function(X,strata,pik,comment=TRUE,method=1)
{
strat=cleanstrata(strata)
H=max(strata)
N=dim(X)[1]
pikstar=rep(0,times=N)
for(h in 1:H) 
{
if(comment==TRUE) cat("\nFLIGHT PHASE OF STRATUM",h)
pikstar[strata==h]=fastflightcube(cbind(X[strata==h,],pik[strata==h]),pik[strata==h],1,comment) 
}
if(comment==TRUE) cat("\nFINAL TREATMENT")
XN=cbind(disjunctive(strata)*pik,X)/pik*pikstar
if(is.null(colnames(X))==FALSE)
    colnames(XN)<-c(paste("Stratum", 1:H, sep = ""),colnames(X))
samplecube(XN,pikstar,1,comment,method) 
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/balancedtwostage.R"
"balancedtwostage" <-
function(X,selection,m,n,PU,comment=TRUE,method=1)
{
N=dim(X)[1]
p=dim(X)[2]
str=cleanstrata(PU)
M=max(PU)
res1=balancedcluster(X,m,PU,method,comment)
if(selection==2) 
        {
         pik2=rep(n/N*M/m,times=N);
         if(n/N*M/m>1) stop("at the second stage, inclusion probabilities larger than 1");
         }
if(selection==1) 
        {
        pik2=inclusionprobastrata(str,rep(n/m ,times=max(str)));
        if(max(pik2)>1) stop("at the second stage, inclusion probabilities larger than 1");
        }
liste=(res1[,1]==1)
sf=rep(0,times=N)
sf[liste]=balancedstratification(array(X[liste,]/res1[,2][liste],c(sum(as.integer(liste)),p)),cleanstrata(str[liste]),pik2[liste],comment,method)
x=cbind(sf,res1[,2]*pik2,res1[,1],res1[,2],pik2)
colnames(x)=c("second_stage","final_pik", "primary","pik_first_stage", "pik_second_stage")
x 
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/calib.r"
calib<- function(Xs, d , total, q=rep(1, length(d)), method=c("linear", "raking", "truncated", "logit"), bounds=c(low=0,upp=10),description=FALSE,max_iter=500)
{
if (any(is.na(Xs)) | any(is.na(d)) | any(is.na(total)) | any(is.na(q)))
      stop("the input should not contain NAs")
if(!(is.matrix(Xs) | is.array(Xs))) Xs=as.matrix(Xs)
if(is.matrix(Xs))
  if(length(total)!=ncol(Xs))
  stop("Xs and total have different dimensions") 
if(is.vector(Xs) & length(total)!=1)
	stop("Xs and total have different dimensions") 
if(any(is.infinite(q))) stop("there are Inf values in the q vector")
if(missing(method)) stop("Specify a method")
if(!(method %in% c("linear", "raking", "logit", "truncated"))) 
stop("the specified method is not in the list")
if(method %in% c("linear", "raking") & !missing(bounds))
stop("for the linear and raking the bounds are not allowed")
require(MASS)
EPS=.Machine$double.eps
EPS1=1e-6
n=length(d)
lambda=as.matrix(rep(0,n))
lambda1=ginv(t(Xs*d*q)%*%Xs,tol=EPS)%*%(total-c(t(d)%*%Xs))
if(method=="linear" | max(abs(lambda1))<EPS) 
          g=1+q*c(Xs%*%lambda1)
else
if(method=="truncated")
      {
     if(!missing(bounds))
       {if(bounds[2]<=1 || bounds[1]>=1 || bounds[1]>bounds[2]) stop("The conditions low<1<upp are not fulfilled")}
     else stop("Specify the bounds")
     g=1+q*c(Xs%*%lambda1)
     list=1:length(g)
     l=0
     g1=g
     Xs1=Xs; d1=d; t2=total;list1=1:length(g);q1=q
     while(l<max_iter)  
        { l=l+1
          if(any(g<bounds[1]) | any(g>bounds[2]))  
               {
		    g[g<bounds[1]]=bounds[1] 
		    g[g>bounds[2]]=bounds[2]
                list=(1:length(g))[g>bounds[1] & g<bounds[2]]
                if(length(list)!=0)
                {
                g1=g[list]
                t2=total-c(t(g[-list]*d[-list])%*%Xs[-list,])
                Xs1=Xs[list,]
                d1=d[list]
                q1=q[list]
                list1=list  
                }
                }
                t1=c(t(d1)%*%Xs1) 
                lambda1=ginv(t(Xs1*d1*q1)%*%Xs1,tol=EPS)%*%(t2-t1)
	          if(length(list1)>1)
         	    g1=1+q1*c(Xs1%*%lambda1)
                else if(length(list1)==1) 
                     {g1=1+q1*c(as.vector(Xs1)%*%as.vector(lambda1))}
                g[list1]=g1
                tr=crossprod(Xs,g*d)
                if(max(abs(tr-total)/total)<EPS1 & all(g>=bounds[1] & g<=bounds[2])) break; 
          
                      
	    } 
  if(l==max_iter) {
                   cat("Not convergence in", max_iter, "iterations with the given bounds. \n");
                   cat("The bounds for the g-weights are:", min(g), " and ", max(g),"\n")   
                   cat(" and the g-weights are given by g\n")
                   }
     } 
else 
if(method=="raking")
{
lambda=as.matrix(rep(0,ncol(Xs)))
w1=as.vector(d*exp(Xs%*%lambda*q))
T=t(Xs)
for(l in 1:max_iter)
	{
      phi=t(Xs)%*%w1-total
	T1=NULL
      for(k in 1:nrow(T)) T1=rbind(T1,T[k,]*w1)
      phiprim=T1%*%Xs 
	lambda=lambda-ginv(phiprim,tol=EPS)%*%phi
     	w1=as.vector(d*exp(Xs%*%lambda*q))
      if(any(is.na(w1))) {warning("Not convergence");g=NULL;break}
      tr=crossprod(Xs,w1)
      if(max(abs(tr-total)/total)<EPS1) break; 
        }
   if(l==max_iter) {warning("Not convergence");g=NULL;}
   else g=w1/d
}
else if(method=="logit")
{
if(bounds[2]<=1 || bounds[1]>=1 || bounds[1]>bounds[2]) stop("The conditions low<1<upp are not fulfilled")
A=(bounds[2]-bounds[1])/((1-bounds[1])*(bounds[2]-1))
u=rep(1,length(d))
F=(bounds[1]*(bounds[2]-1)+bounds[2]*(1-bounds[1])*u)/(bounds[2]-1+(1-bounds[1])*u)
w1=as.vector(d*F)
T=t(Xs)
for(k in 1:nrow(T)) T[k,]=T[k,]*w1
phiprim=ginv(T%*%Xs,tol=EPS) 
g=F  
tr=crossprod(Xs, w1)
if(max(abs(tr-total)/total)>EPS1 | any(g<bounds[1]) | any(g>bounds[2]))  
{
lambda1=rep(0,ncol(Xs))
list=1:length(g)
t2=total
Xs1=Xs
d1=d 
#D1=D
g1=g
q1=q
list1=1:length(g)
for(l in 1:max_iter)
{
if(any(g<bounds[1]) | any(g>bounds[2]))  
	{g[g<bounds[1]]=bounds[1]
       g[g>bounds[2]]=bounds[2] 
       list=(1:length(g))[g>bounds[1] & g<bounds[2]]    
       if(length(list)!=0)
                {
                g1=g[list]
                t2=total-c(t(g[-list]*d[-list])%*%Xs[-list,])
                Xs1=Xs[list,]
                d1=d[list]
                q1=q[list]
                list1=list  
                }
       }
t1=c(t(d1)%*%Xs1) 
phi=t(Xs1)%*%as.vector(d1*g1)-t1
T=t(Xs1)
for(k in 1:nrow(T))  T[k,]=T[k,]*d1*g1 
phiprime=T%*%Xs1
lambda1=lambda1-ginv(phiprime,tol=EPS)%*%(as.vector(phi)-t2+t1)
u=exp(A*(Xs1%*%lambda1*q1))
F=g1=(bounds[1]*(bounds[2]-1)+bounds[2]*(1-bounds[1])*u)/(bounds[2]-1+(1-bounds[1])*u)  
if(any(is.na(g1))) {warning("not convergence");g1=g=NULL;break}
g[list1]=g1
tr=crossprod(Xs,g*d)
if(max(abs(tr-total)/total)<EPS1 & all(g>=bounds[1] & g<=bounds[2])) break; 
}
if(l==max_iter) {
                 cat("not convergence in", max_iter, "iterations with the given bounds. \n");
                 cat("the bounds for the g-weights are:", min(g), " and ", max(g),"\n")   
                 cat(" and the g-weights are given by g\n")
                 g=NULL
                     }
}
}
if(description && !is.null(g)) {
par(mfrow=c(3,2),pty="s")
ratio=1/g 
hist(ratio);boxplot(ratio,main="Boxplot of ratio w/d")
hist(d);boxplot(d,main="Boxplot of d")
hist(g*d);boxplot(g*d,main="Boxplot of w=g*d")
if(method %in% c("truncated","raking","logit"))
cat("number of iterations ",l,"\n")
cat("summary - initial weigths d\n")
print(summary(d))
cat("summary - final weigths w=g*d\n")
print(summary(as.vector(g*d)))
}
g
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/calibev.r"
calibev<-function(Ys,Xs,total,pikl,d,g,q=rep(1,length(d)),with=FALSE,EPS=1e-6)
{
if(any(is.na(g))) 
   stop("There are missing values in g")
stopifnot((ns <- length(g)) >= 1)
if(min(pikl)==0) 
{ss=NULL
warning("There are zero values in the 'pikl' matrix. The variance estimator can not be computed.\n")
}
piks=as.vector(diag(pikl))
if(!checkcalibration(Xs,d,total,g,EPS)$result) 
stop("The calibration is not possible. The calibration estimator is not computed.\n")
if(is.data.frame(Xs)) Xs=as.matrix(Xs)
if(!is.vector(Ys)) Ys=as.vector(Ys)
if(is.matrix(Xs)) n=nrow(Xs)
else n=length(Xs)
if(ns!=length(Ys) | ns!=length(piks) | ns!=n | ns!=length(d)) stop("The parameters have different sizes.\n")
w=g*d
wtilde=w*q
B=t(Xs*wtilde) #this is for t(Xs)%*%diag(wtilde)
beta=ginv(B%*%Xs)%*%B%*%Ys
e=Ys-Xs%*%beta
ss<-0
for(k in 1:ns)
 {ss2<-0
 for(l in 1:ns)
       if(!with) 
       ss2<-ss2+(1-piks[k]*piks[l]/pikl[k,l])*w[k]*e[k]*w[l]*e[l]
       else 
       ss2<-ss2+(1-piks[k]*piks[l]/pikl[k,l])*d[k]*e[k]*d[l]*e[l]
 ss<-ss+ss2
 }
list(calest=sum(w*Ys),evar=as.numeric(ss))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/checkcalibration.R"
checkcalibration<-function (Xs, d, total, g, EPS=1e-6) 
{
 if(is.null(g)) stop("the g-weight vector is null")
 if(!is.matrix(Xs)) Xs<-as.matrix(Xs)
 tr <- crossprod(Xs, g*d)
 if(max(abs(tr-total)/total)<EPS) {result<-TRUE
                           message<-"the calibration is done"
                           value<-EPS  
                           }
else {message<-cat("the calibration cannot be done. The max EPS value is given by 'value'.\n")
      value<-max(abs(tr-total)/total)
      result<-FALSE}
list(message=message,result=result,value=value)
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/cleanstrata.R"
"cleanstrata" <-
function(strata)
{
a=sort(unique(strata)) 
b=1:length(a) 
names(b)=a 
as.vector(b[as.character(strata)]) 
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/cluster.r"
cluster<-function(data, clustername, size, method=c("srswor","srswr","poisson","systematic"),pik,description=FALSE)
{
if(size==0) stop("the size is zero")
if(missing(method)) {warning("the method is not specified; by default, the method is srswor")
                     method="srswor"
                    }
if(!(method %in% c("srswor","srswr","poisson","systematic"))) 
  stop("the name of the method is wrong")
if(method %in% c("poisson","systematic") & missing(pik)) stop("the vector of probabilities is missing")
data=data.frame(data)
index=1:nrow(data)
if(missing(clustername))
   {
   if(method=="srswor")
	result=data.frame(index[srswor(size,nrow(data))==1],rep(size/nrow(data),size))
   if(method=="srswr")
     {
      s=srswr(size,nrow(data))
      st=s[s!=0]
      l=length(st)
      result=data.frame(index[s!=0])
      if(size<=nrow(data))    
          result=cbind.data.frame(result,st,prob=rep(size/nrow(data),l))
      else 
          {prob=rep(size/nrow(data),l)/sum(rep(size/nrow(data),l))
           result=cbind.data.frame(result,st,prob)
          }  
  
      colnames(result)=c("ID_unit","Replicates","Prob")
     }
   if(method=="poisson")
    {
     pikk=inclusionprobabilities(pik,size)
     s=(UPpoisson(pikk)==1)
     if(length(s)>0)
     result=data.frame(index[s],pikk[s])
     if(description) 
      cat("\nPopulation total and number of selected units:",nrow(data),length(s),"\n") 
     }
    if(method=="systematic")
    {
     pikk=inclusionprobabilities(pik,size)
     s=(UPsystematic(pikk)==1)
     result=data.frame(index[s],pikk[s])
     }
   if(method!="srswr")
          colnames(result)=c("ID_unit","Prob")
   if(description) cat("\nPopulation total and number of selected units:",nrow(data),sum(size),"\n")
   }
else
{
data=data.frame(data)
m=match(clustername,colnames(data))
if(length(m)>1) stop("there are too many specified variables as clusters")
if(is.na(m)) stop("the cluster name is wrong")
x1=as.factor(data[,m])
result=NULL
if(nlevels(x1)==0) stop("the cluster name has 0 modalities")
else
{
nr_cluster=nlevels(x1)
if(method=="srswor")
 {
s=as.data.frame(levels(x1)[srswor(size,nr_cluster)==1])
names(s)=c("cluster")
r=cbind.data.frame(index,data[,m])
names(r)=c("index","cluster")
r=merge(r,s,by.x="cluster",by.y="cluster",sort=TRUE)
result=cbind.data.frame(r,rep(size/nr_cluster,nrow(r)))
 }
if(method=="srswr")
{
s=srswr(size,nr_cluster)
st=cbind.data.frame(levels(x1)[s!=0],s[s!=0])
names(st)=c("cluster","repl")
r=cbind.data.frame(index,data[,m])
names(r)=c("index","cluster")
r=merge(r,st,by.x="cluster",by.y="cluster") 
if(size<=nr_cluster)  
 result=cbind.data.frame(r,rep(size/nr_cluster,nrow(r)))
else
{prob=rep(size/nr_cluster,nrow(r))/sum(rep(size/nr_cluster,nrow(r))) 
result=cbind.data.frame(r,prob)
}
}
if(method=="systematic")
{pikk=inclusionprobabilities(pik,size)
 s=(UPsystematic(pikk)==1)  
 st=cbind.data.frame(levels(x1)[s],pikk[s])
 names(st)=c("cluster","prob")
 r=cbind.data.frame(index,data[,m])
 names(r)=c("index","cluster")
 result=merge(r,st,by.x="cluster",by.y="cluster") 
}
if(method=="poisson") 
{pikk=inclusionprobabilities(pik,size)
 s=(UPpoisson(pikk)==1)
 if(any(s))
             { 
 st=cbind.data.frame(levels(x1)[s],pikk[s])
 names(st)=c("cluster","prob")
 r=cbind.data.frame(index,data[,m])
 names(r)=c("index","cluster")
 result=merge(r,st,by.x="cluster",by.y="cluster") 
 if(description) 
  {cat("Number of selected clusters:" ,sum(s),"\n")
   cat("\nPopulation total and number of selected units:",nrow(data),nrow(result),"\n")
  }
               }
             else 
		{if(description) 
		 {cat("Number of selected clusters: 0\n")
		  cat("Population total and number of selected units:",nrow(data),0,"\n")
              }
             result=NULL 
     	       }
}
if(method=="srswr")
{
colnames(result)=c(clustername,"ID_unit","Replicates","Prob")
if(description)
  {cat("Number of selected clusters:",length(s[s!=0]),"\n")
   cat("Population total and number of selected units", nrow(data),nrow(result),"\n")
   }
}
else
if(!is.null(result))
colnames(result)=c(clustername,"ID_unit","Prob")
if(description & !(method %in%c("poisson","srswr")))  
      {cat("Number of selected clusters:",size,"\n")
       cat("Population total and number of selected units", nrow(data),nrow(result),"\n")
       }
}
}
result
}




#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/disjunctive.R"
"disjunctive" <-
function(strata)
{ss=cleanstrata(strata)
m=matrix(0,length(strata),length(unique(strata)))
for(i in 1:length(ss)) m[i,ss[i]]=1
m
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/fastflightcube.R"
"fastflightcube" <-
function(X,pik,order=1,comment=TRUE) 
{ 
EPS = 1e-11
"algofastflightcube" <-
function(X,pik) 
{ 

"jump" <-
function(X,pik){ 
N = length(pik)
p = round(length(X)/length(pik))
X<-array(X,c(N,p))
X1=cbind(X,rep(0,times=N)) 
kern<-svd(X1)$u[,p+1] 
listek=abs(kern)>EPS
buff1<-(1-pik[listek])/kern[listek]
buff2<- -pik[listek]/kern[listek]
la1<-min( c(buff1[(buff1>0)] , buff2[(buff2>0)]) ) 
pik1<- pik+la1*kern 
buff1<- -(1-pik[listek])/kern[listek] 
buff2<- pik[listek]/kern[listek]
la2<-min(c(buff1[(buff1>0)] , buff2[(buff2>0)])) 
pik2<- pik-la2*kern 
q<-la2/(la1+la2)  
if (runif(1)<q) pikn<-pik1 else pikn<-pik2 
pikn 
}

N = length(pik)
p = round(length(X)/length(pik))
X<-array(X,c(N,p));
A<- X/pik; 
B<-A[1:(p+1),]; 
psik <- pik[1:(p+1)]; 
ind<-seq(1,p+1,1); 
pp=p+2; 
B<-array(B,c(p+1,p)); 
while(pp<=N) 
{ 
psik <- jump(B,psik); 
liste<- (psik>(1-EPS) | psik<EPS); 
i<- 0; 
while(i <=(p) & pp<=N) 
{ i=i+1;if(liste[i]==TRUE) 
{pik[ind[i]]=psik[i]; 
psik[i]=pik[pp] ; 
B[i,]=A[pp,];B<-array(B,c(p+1,p)); 
ind[i]=pp; 
pp=pp+1 } 
} 
} 
if(length(pik[(pik>EPS & pik<(1-EPS))])==(p+1)) psik <- jump(B,psik); 
pik[ind]=psik; 
pik 
}

"reduc" <-
function(X)
{
EPS=0.0000000001
N=dim(X)[1]
Re=svd(X)
array(Re$u[,(Re$d>EPS)] , c(N,sum(as.integer(Re$d>EPS))))
}

N = length(pik);
p = round(length(X)/length(pik));
X<-array(X,c(N,p));
if (order==1) o<-sample(N,N)  else
   {
   if(order==2) o<-seq(1,N,1) 
    else o<-order(pik,decreasing=TRUE)
    }
liste<-o[(pik[o]>EPS & pik[o]<(1-EPS))];
if(comment==TRUE){
cat("\nBEGINNING OF THE FLIGHT PHASE\n");
cat("The matrix of balanced variable has",p," variables and ",N," units\n");
cat("The size of the inclusion probability vector is ",length(pik),"\n");
cat("The sum of the inclusion probability vector is ",sum(pik),"\n");
cat("The inclusion probability vector has ",length(liste)," non-integer elements\n");
} 
pikbon<-pik[liste]; 
Nbon=length(pikbon); 
Xbon<-array(X[liste,] ,c(Nbon,p)); 
pikstar<-pik; 
flag=0;
if(Nbon>p){if(comment==TRUE) cat("Step 1  ");
           pikstarbon<-algofastflightcube(Xbon,pikbon); 
           pikstar[liste]=pikstarbon; 
           flag=1
           }
liste<-o[(pikstar[o]>EPS & pikstar[o]<(1-EPS))];
pikbon<-pikstar[liste]; 
Nbon=length(pikbon); 
Xbon<-array(X[liste,] ,c(Nbon,p)); 
pbon=dim(Xbon)[2]
if(Nbon>0){
          Xbon=reduc(Xbon)
          pbon=dim(Xbon)[2]
          }
k=2
while(Nbon>pbon & Nbon>0){
           if(comment==TRUE) cat("Step ",k,",  ");
           k=k+1;
           pikstarbon<-algofastflightcube(Xbon/pik[liste]*pikbon,pikbon); 
           pikstar[liste]=pikstarbon; 
           liste<-o[(pikstar[o]>EPS & pikstar[o]<(1-EPS))];
           pikbon<-pikstar[liste]; 
           Nbon=length(pikbon); 
           Xbon<-array(X[liste,] ,c(Nbon,p)); 
           if(Nbon>0)
               {
               Xbon=reduc(Xbon)
               pbon=dim(Xbon)[2]
               }
           flag=1
           }
if(comment==TRUE) if(flag==0) cat("NO FLIGHT PHASE");
if(comment==TRUE) cat("\n")
pikstar 
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/getdata.r"
getdata<-function(data,m)
{
if(!is.data.frame(data)) data=as.data.frame(data)
if(is.null(names(data))) stop("the column names are missing")
if(is.vector(m)& !is.list(m))
   {res=NULL 
    if(is.null(names(m)))
      if(all(m %in% c(0,1))) {
      res=NULL
   if(!("ID_unit" %in% names(data))) 
   {res=cbind.data.frame((1:length(m))[m==1],data[m==1,])
   names(res)=c("ID_unit",names(data))
   }
   else res=data[m==1,]
                                 }
      else  
   res=data[m,]
    }
else
if(is.data.frame(m))
     {res=NULL 
      if(!is.null(names(m)))
		{   
      mm=match(names(data),names(m),nomatch=0)
      index=(1:ncol(data))[mm==0]
      if(length(index)>0)
	      {res=cbind.data.frame(data[m$"ID_unit",index],m)
             names(res)[1:length(index)]=names(data)[index] 
 		}
          }
     }
else
if(is.list(m))
 {res=list()
  if(length(m)>=1)
  for(j in 1:length(m))
	{
     mm=match(names(data),names(m[[j]]),nomatch=0)
     index=(1:ncol(data))[mm==0]
     if(length(index)>0)
     {
     res[[j]]=cbind.data.frame(data[m[[j]]$"ID_unit",index],m[[j]])
     names(res[[j]])[1:length(index)]=names(data)[index] 
      }}
  else res=m
} 
res
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/inclusionprobabilities.R"
"inclusionprobabilities" <- function(a,n)
{   nnull = length(a[a == 0])
    nneg = length(a[a < 0])
    if (nnull > 0) 
        warning("there are zero values in the initial vector a\n")
    if (nneg > 0) {
        warning("there are ", nneg, " negative value(s) shifted to zero\n")
        a[(a < 0)] = 0
    }
    if(identical(a,rep(0,length(a)))) pik1=a
    else
    {
    pik1 =n * a/sum(a)
    pik=pik1[pik1>0]
    list1=pik1>0
    list = pik >= 1
    l = length(list[list == TRUE])
    if(l>0)
    {
    l1=0
    while (l != l1) {
	x=pik[!list]
      x=x/sum(x)
      pik[!list] = (n - l)*x
      pik[list] = 1
      l1 = l
      list = (pik >= 1)
      l = length(list[list == TRUE])
    }
    pik1[list1]=pik
    }
    }
pik1
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/inclusionprobastrata.R"
"inclusionprobastrata" <-
function(strata,nh)
{
EPS=1e-6
N=length(strata)
v=unique(strata)
ma=max(v)
if(min(v)<1) stop("the stratification variable has incorect values (less than 1)\n")
Nh=rep(0,times=ma)
for(i in 1:N) Nh[strata[i]]=Nh[strata[i]]+1
for(i in 1:ma) 
    if(nh[i]/Nh[i]>1+EPS) 
        warning("in stratum ",i," the sample size is larger than the population size\n")
pik=rep(0,times=max(v))
for(i in 1:N)  pik[i]=nh[strata[i]]/Nh[strata[i]]
pik
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/landingcube.R"
"landingcube" <-
function(X,pikstar,pik,comment=TRUE) 
# landing phase of the cube method
######################################################
{ 
# extraction of the non-integer values for the landing phase 
require(MASS)
require(lpSolve)
EPS=1e-11
p=dim(X)[2]
N=dim(X)[1]
liste=(pikstar>EPS & pikstar<(1-EPS)) 
pikland=pikstar[liste] 
Nland=length(pikland) 
Xland=array(X[liste,] ,c(Nland,p)) 
nland=sum(pikland) 
FLAGI=(abs(nland-round(nland))<EPS);
# construction of the list of possible samples 
if(comment==TRUE){
cat("\n\nBEGINNING OF THE LANDING PHASE\n");
cat("At the end of the flight phase, there remain ",Nland,"non integer probabilities","\n");
cat("The sum of these probabilities is ",nland,"\n");
cat("This sum is "); 
if(FLAGI)  cat(" integer\n") else cat(" non-integer\n");
}
if(FLAGI) 
     {
      pikland=round(nland)*pikland/nland;
      nland=round(nland);
      SSS= writesample(nland,Nland) 
      }
else 
      SSS= rbind( writesample(trunc(nland),Nland), 
           writesample(trunc(nland)+1,Nland) ); 
lll=nrow(SSS);
if(comment==TRUE){
cat("The linear program will consider ",lll," possible samples\n")
} 
# computation of the cost 
Asmp=matrix(0,p,lll);
for(i in 1:lll)  {Asmp[,i]=t(Xland/pik[liste]) %*% (SSS[i,]-pikland) }
A<- X[pik>EPS,]/pik[pik>EPS]; 
cost=rep(0,times=lll) 
for(i in 1:lll) 
cost[i]=t(Asmp[,i]) %*% ginv(t(A) %*% A) %*% Asmp[,i] 
# linear programming
V = t(cbind(SSS,rep(1,times=lll)))
b=c(pikland,1)
constdir=rep("==",times=(Nland+1))
x=lp("min",cost,V,constdir,b)$solution
# choice of the sample
u=runif(1,0,1)
i=0
ccc=0
while(ccc<u) {i=i+1;ccc=ccc+x[i]}
if(comment==TRUE)
    {
     cat("The mean cost is ",mean(cost),"\n")
     cat("The smallest cost is ",min(cost),"\n")
     cat("The largest cost is ",max(cost),"\n")
     cat("The cost of the selected sample is",cost[i])
     }
pikfin=pikstar
pikfin[liste]=SSS[i,];
pikfin
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/mstage.r"
mstage<-function(data, stage=c("stratified","cluster",""),varnames, size, method=c("srswor","srswr","poisson","systematic"), pik, description=FALSE)
{
if(missing(size)) stop("the size argument is missing")
if(!missing(stage) & missing(varnames)) stop("please indicate the stage argument")
if(!missing(stage)) 
  {number=length(stage)
   for(i in 1:length(stage))
    if(!(stage[i] %in% c("stratified","cluster",""))) stop("the stage argument is wrong")
  }
else number=length(size)
if(number>1) 
	{if(!missing(varnames))
	 {  
       if(!is.list(size)) stop("the size must be a list")
       size=as.list(size)
       varnames=as.list(varnames)
	 size1=size[[1]] 
       varnames1=varnames[[1]] 
       if(method %in% c("systematic","poisson"))
	       pik1=pik[[1]]
	 }
	else {size1=size[[1]];varnames1=NULL
           if(method %in% c("systematic","poisson")) pik1=pik[[1]]}
      }
else {size1=size
      if(method %in% c("systematic","poisson")) pik1=pik}
if(description)
 cat("STAGE 1","\n")
if(missing(stage))
 {if(missing(varnames))
  if(missing(method))
 	s=strata(data, stratanames=NULL, size1, description)
  else 
      if(method %in% c("systematic","poisson")) 
 	s=strata(data, stratanames=NULL, size=size1, method, pik1, description)
      else s=strata(data, stratanames=NULL, size=size1, method, description)
  else s=strata(data, stratanames=NULL, size1, method, pik1, description)
  }
else 
 if(stage[1]=="stratified") 
  {s=strata(data,varnames1,size1, method, pik=pik1, description)
   dimension_st=table(s$Stratum)
   }
 else {
       s=cluster(data, varnames1, size1, method, pik1, description)
       if(is.null(s)) stop("0 selected units in the first stage")
       m=match(varnames1,names(s))
       nl=nlevels(as.factor(s[,m]))
       lev=levels(as.factor(s[,m]))
       if(nl>=1)
       {
       dimension_cl=numeric(nl)
       for(i in 1:nl)
         dimension_cl[i]=nrow(subset(s,as.factor(s[,m])==lev[i]))
       } 
       dimension=dimension_cl
       }    
if(is.null(s)) stop("0 selected units in the first stage")
if(number>1)
 if((is.element("cluster",stage) & is.element("stratified",stage)))
  result=getdata(data,s)
  else result=s
res=list()
res[[1]]=s

if(number>=2)
for(j in 2:number)
{
if(description) cat("STAGE ",j,"\n")
if(!missing(varnames))
{ 
if(stage[[j]]=="cluster")
{
if(stage[[j-1]]=="stratified") 
{
k=length(dimension_st)            
s1=NULL  
limit=0
dimension=list()
if(k>=1)
for(ii in 1:k)
                {
                 r=res[[j-1]][(limit+1):(limit+dimension_st[ii]),]
                 r=getdata(data,r) 
                 m=match(varnames[[j]],names(r))
                 if(method %in% c("systematic","poisson")) 
                 {index=res[[j-1]][(limit+1):(limit+dimension_st[ii]),]$ID_unit 
                 pikk=pik[[j]][index]
                 if(!is.null(r))
		     s3=cluster(r,clustername=varnames[[j]],size=size[[j]][ii],method=method, pik=pikk,description)                        
                 else s3=NULL 
                 }
                 else s3=cluster(r,clustername=varnames[[j]],size=size[[j]][ii],method=method, pik,description)                        
                limit=limit+dimension_st[ii]          
                if(method=="srswr")
                {s3=cbind.data.frame(r[s3$ID_unit,m],r[s3$ID_unit,]$ID_unit,s3$Replicates,s3$Prob,r[s3$ID_unit,]$Prob*s3$Prob)
                 colnames(s3)=c(varnames[[j]],"ID_unit","Replicates",paste("Prob_",j,"_stage"),"Prob")
                 }
		   else if(!is.null(s3))
		    {s3=cbind.data.frame(r[s3$ID_unit,m],r[s3$ID_unit,]$ID_unit,s3$Prob,result[s3$ID_unit,]$Prob*s3$Prob)
		     colnames(s3)=c(varnames[[j]],"ID_unit",paste("Prob_",j,"_stage"),"Prob")
		    }
if(!is.null(s3))
{
m=match(varnames[[j]],names(s3))
for(l in 1:nlevels(as.factor(s3[,m])))
	dimension=c(dimension,table(s3[,m])[l])
s1=rbind(s1,s3)
}
}
}
else
 if(stage[[j-1]]=="cluster") 
{
k=length(dimension_cl)            
s1=NULL  
limit=0
dimension=list()
if(k>length(size[[j]])) 
           {warning("the number of selected clusters at the previous stage is larger than the size argument")                  
            warning("the size 1 is added")
            size1=size[[j]]
            for(i in 1:(k-length(size[[j]])))
            size1=c(size1,1)
            }
else size1=size[[j]]
if(k>=1)
for(ii in 1:k)
                {
                 r=res[[j-1]][(limit+1):(limit+dimension_cl[ii]),]
                 r=getdata(data,r) 
                 m=match(varnames[[j]],names(r))
                 if(method %in% c("systematic","poisson")) 
                 { 
			m1=match(varnames[[j-1]],names(r))
			m2=match(varnames[[j-1]],names(data))
			mm=match(r[1,m1],levels(factor(data[,m2])))
                  pikk=as.numeric(pik[[j]][[mm]])  
                  if(!is.null(r)) 
		      s3=cluster(r,clustername=varnames[[j]],size=size1[[ii]],method=method, pik=pikk,description)                        
                  else s3=NULL
                  }
                 else s3=cluster(r,clustername=varnames[[j]],size=size1[ii],method=method, pik,description)                        
                limit=limit+dimension_cl[ii]           
                if(method=="srswr")
                {s3=cbind.data.frame(r[s3$ID_unit,m],r[s3$ID_unit,]$ID_unit,s3$Replicates,s3$Prob,r[s3$ID_unit,]$Prob*s3$Prob)
                 colnames(s3)=c(varnames[[j]],"ID_unit","Replicates",paste("Prob_",j,"_stage"),"Prob")
                 }
		   else if(!is.null(s3))
		    {s3=cbind.data.frame(r[s3$ID_unit,m],r[s3$ID_unit,]$ID_unit,s3$Prob,result[s3$ID_unit,]$Prob*s3$Prob)
		     colnames(s3)=c(varnames[[j]],"ID_unit",paste("Prob_",j,"_stage"),"Prob")
		    }
if(!is.null(s3))
{
m=match(varnames[[j]],names(s3))
for(l in 1:nlevels(as.factor(s3[,m])))
	dimension=c(dimension,table(s3[,m])[l])
s1=rbind(s1,s3)
}
} 
}
} 
else       
		  if(j>1)
               {k=length(dimension) 
                s1=NULL  
                limit=0
                count=0 
                if(k>length(size[[j]])) 
                  {warning("the number of selected clusters at the previous stage is larger than the size argument")                  
                   warning("the size 1 is added")
                   size1=size[[j]]
                   for(i in 1:(k-length(size[[j]])))
                   	size1=c(size1,1)
                   }
                else size1=size[[j]]
                if(k>=1)
                for(i in 1:k)
                 for(ii in 1:length(dimension[[i]]))
                {
                 r=res[[j-1]][(limit+1):(limit+dimension[[i]][ii]),]
                 count=count+1
                 if(method %in% c("systematic","poisson")) 
                 {index=res[[j-1]][(limit+1):(limit+dimension[[i]][ii]),]$ID_unit 
                 pikk=pik[[j]][index]
                 if(!is.null(r))
                 s2=strata(r,NULL,size=size1[count],method=method, pik=pikk,description)                        
                 else s2=NULL
                 }
                else s2=strata(r,NULL,size=size1[count],method=method, pik,description)                        
                limit=limit+dimension[[i]][ii]              
                if(method=="srswr")
                {s2=cbind.data.frame(r[s2$ID_unit,]$ID_unit,s2$Replicates,s2$Prob,r[s2$ID_unit,]$Prob*s2$Prob)
                 colnames(s2)=c("ID_unit","Replicates",paste("Prob_",j,"_stage"),"Prob")
                 }
 		  else if(!is.null(s2))
		  {
		   s2=cbind.data.frame(r[s2$ID_unit,]$ID_unit,s2$Prob,r[s2$ID_unit,]$Prob*s2$Prob)
		   colnames(s2)=c("ID_unit",paste("Prob_",j,"_stage"),"Prob")
   
		   }
if(!is.null(s2))
s1=rbind(s1,s2)
                 } 
}
}		  
else 
  {
 if(missing(stage))
 {if(missing(method))
 	s1=strata(result, stratanames=NULL, size=size[[j]], description=description)
  else  
if(method=="poisson" | method=="systematic") 
   s1=strata(result,stratanames=NULL, size=size[[j]], method=method, pik=pik[[j]], description=description)
else 
  s1=strata(result,stratanames=NULL,size=size[[j]], method=method, description=description)
  if(method=="srswr")
  {s1=cbind.data.frame(result[s1$ID_unit,]$ID_unit,s1$Replicates,s1$Prob,result[s1$ID_unit,]$Prob*s1$Prob)
   colnames(s1)=c("ID_unit","Replicates",paste("Prob_",j,"_stage"),"Prob")
  }
  else
  {
   s1=cbind.data.frame(result[s1$ID_unit,]$ID_unit,s1$Prob,result[s1$ID_unit,]$Prob*s1$Prob)
   colnames(s1)=c("ID_unit",paste("Prob_",j,"_stage"),"Prob")
   }
  }
}
if(!is.null(s1)) 
	{result=s1
	res[[j]]=result}
else number=number-1
}
if(!is.null(names(res[[1]])))
  {m=match("Prob", names(res[[1]]))    
   names(res[[1]])[m]="Prob_ 1 _stage" 
   }
names(res)=c(1:number)
res
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/postest.r"
postest<-function(data, y, pik, NG, description=FALSE)
{
if (missing(data) | missing(y) | missing(pik) | missing(NG)) 
    stop("incomplete input")
str <- function(st, h, n) .C("str", as.double(st), as.integer(h), 
        as.integer(n), s = double(n), PACKAGE = "sampling")$s
data=as.data.frame(data)
sample.size=nrow(data)
t_post=0
if(!is.null(colnames(data)))
{m = match("Stratum", colnames(data))
if(!any(is.na(m))) 
        {
          m = match("poststratum", colnames(data))
          if (any(is.na(m))) 
            stop("the column 'poststratum' is missing")
        h=unique(data$Stratum)
        g=unique(data$poststratum)
        for(j in 1:length(g))   
         {p=str(data$poststratum, g[j], sample.size)
          Ng=sum(NG[,j])
          t1=t2=0
          for (i in 1:length(h))         
            {s = str(data$Stratum, h[i], sample.size)
             shg=s*p
             if(!all(shg==0))  
                {
                nhg=length(shg[shg==1])
                t1=t1+sum(y[shg==1]/pik[shg==1])
                t2=t2+sum(1/pik[shg==1])               
                if(description)
                     {cat("Stratum ",j,", postratum ", i," \n")
                      cat("the postratified estimator is:",Ng*t1/t2,"\n")
                      }
                t_post=t_post+Ng*t1/t2
               }
             else if(description) 
                cat("Stratum ",j,", postratum ", i," empty intersection set \n")
             }}
    }
else {
      g=unique(data$poststratum)
      for(j in 1:length(g))   
         {p=str(data$poststratum, g[j], sample.size)
          Ng=NG[j]
          t1=Ng*sum(y[p==1]/pik[p==1])/sum(1/pik[p==1])
          t_post=t_post+t1
          if(description)
                     {cat("postratum ", j," \n")
                      cat("the postratified estimator is:",t1,"\n")
                     }
          }
       }
}
else  stop("the column names in data are missing")  
t_post
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/poststrata.r"
poststrata<-function(data, postnames = NULL) 
{
 if (missing(data) | missing(postnames)) stop("incomplete input")
 data = data.frame(data)
 if(is.null(colnames(data))) stop("the column names in data are missing")  
 index = 1:nrow(data)
 m = match(postnames, colnames(data))
 if (any(is.na(m))) 
            stop("the names of the poststrata are wrong")
 data2 = cbind.data.frame(data[, m])
 x1 = data.frame(unique(data[, m]))
 colnames(x1) = postnames
 nr_post=0
 post=numeric(nrow(data))
 nh=numeric(nrow(x1))
 for(i in 1:nrow(x1)) 
    { expr=rep(FALSE, nrow(data2))
      for(j in 1:nrow(data2)) expr[j]=all(data2[j, ]==x1[i, ])
      y=index[expr]
      if(is.matrix(y)) 
        nh[i]=nrow(y)
      else nh[i]=length(y)
      post[expr]=i
      }
result=cbind.data.frame(data,post)
names(result)=c(names(data),"poststratum")
list(data=result, npost=nrow(x1))
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/ratioest.r"
ratioest<-function(y,x,Tx,pik)
{if (any(is.na(pik))) 
        stop("there are missing values in pik")
    if (any(is.na(y))) 
        stop("there are missing values in y")
    if (any(is.na(x))) 
        stop("there are missing values in x")
    if (length(y) != length(pik) | length(x)!=length(pik) | length(x)!=length(y)) 
        stop("y, x and pik have different lengths")
    sum(y/pik)*Tx/sum(x/pik)
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/ratioest_strata.r"
ratioest_strata<-function(y,x,TX_strata,pik,strata,description=FALSE)
{
if (missing(x) | missing(y) | missing(pik) | missing(strata)) stop("incomplete input")
if(!is.vector(x)) x=as.vector(x)
if(length(y)!=length(x)) stop("x and y have different sizes")
str <- function(st, h, n) .C("str", as.double(st), as.integer(h), as.integer(n), s = double(n), PACKAGE = "sampling")$s
sample.size = length(y)
h = unique(strata)
s1 = 0
for (i in 1:length(h)) {
s=str(strata, h[i], sample.size)
ys=y[s==1]
xs=x[s==1]
r=ratioest(ys,xs,TX_strata[h[i]],pik[s==1])
s1 = s1 + r
if(description)
 {cat("Stratum ",h[i],", the ratio estimator is:",r,"\n")
  cat("Number of units:",sum(s),"\n")
    }  
  }
if(description) 
  cat("The ratio estimator is:\n")
s1
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/regest.r"
regest<-function(formula,Tx,weights,pikl,n,sigma=rep(1,length(weights)))
{  
cl <- match.call() 
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "weights"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
pik<-1/w
if(!identical(sigma,rep(1,length(pik)))) w<-w/sigma^2
x <- model.matrix(mt, mf, contrasts)
if(ncol(x)==1) x=as.vector(x)
if (any(is.na(pik))) 
        stop("there are missing values in pik")
if (any(is.na(y))) 
        stop("there are missing values in y")
if (any(is.na(x))) 
        stop("there are missing values in x")
if(is.vector(x))
    {if (length(y) != length(pik) | length(x)!=length(pik) | length(x)!=length(y)) 
        stop("y, x and pik have different lengths")
    }
    else
    if(is.matrix(x)) 
    {if (length(y) != length(pik) | nrow(x)!=length(pik) | nrow(x)!=length(y)) 
      stop("y, x and pik have different sizes")
     if(ncol(x)>2 & length(Tx)!=ncol(x)-1)
             stop("x and Tx have different sizes") 
     }

model<-lm(y~x-1,weights=w)
e<-model$residuals
beta<-model$coefficient
# variance of beta, Sarndal p. 194
delta<-matrix(0,nrow(pikl),ncol(pikl))
for(k in 1:(nrow(delta)-1))
 {for(l in (k+1):ncol(delta))
    delta[l,k]<-delta[k,l]<-1-pikl[k,k]*pikl[l,l]/pikl[k,l] 
 delta[k,k]<-1-pikl[k,k]
 }
delta[nrow(delta),ncol(delta)]<-1-pikl[nrow(delta),ncol(delta)]
j_start<-1
if(is.matrix(x))
    {
    if(all(x[,1]==rep(1,nrow(x)))) j_start<-2 
    xx<-as.matrix(x[,j_start:ncol(x)])
    s<-0 
    for(i in 1:ncol(xx))
          if(j_start==2)
    	    s<-s+sum(beta[i+1]*(Tx[i]-HTestimator(xx[,i],pik))) 
          else s<-s+sum(beta[i]*(Tx[i]-HTestimator(xx[,i],pik))) 
    est<-HTestimator(y,pik)+s
    }
else
est<-HTestimator(y,pik)+sum(beta*(Tx-HTestimator(x,pik)))
V<-t(x*w*e)%*%delta%*%(x*e*w)
inv<-ginv(t(x * w) %*% x)
var_beta<-inv%*%V%*%inv
z<-list()
class(z) <- c("regest")
z$call <- cl
z$formula <- formula
z$x <- x
z$y <- y
z$weights<-w
z$regest<-as.numeric(est)
z$coefficients<-beta
z$std_error<-sqrt(diag(var_beta))
z$t_value<-beta/sqrt(diag(var_beta))
# number of degrees of freedom is number of obs-1 if intercept, and number of obs otherwise 
if(j_start==1)
	z$p_value<-2*(1-pt(z$t_value,n-1))
else z$p_value<-2*(1-pt(z$t_value,n))
z$cov_matrix<-var_beta
z
}


#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/regest_strata.r"
regest_strata<-function(formula,weights,Tx_strata,strata,pikl,sigma=rep(1,length(weights)),description=FALSE)
{
cl <- match.call() 
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "weights"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
x <- model.matrix(mt, mf, contrasts)
str <- function(st, h, n) .C("str", as.double(st), as.integer(h), as.integer(n), s = double(n), PACKAGE = "sampling")$s
sample.size = length(y)
h = unique(strata)
s1 = 0
for (i in 1:length(h)) {
s=str(strata, h[i], sample.size)
ys=y[s==1]
xs=x[s==1,]
r=regest(ys~xs-1,Tx=Tx_strata[h[i]],weights=weights[s==1],pikl=pikl[s==1,s==1],n=length(s[s==1]),sigma[s==1])
est=r$regest
s1 = s1 + est
if(description)
 {cat("Stratum ",h[i],", the regression estimator is:",est,"\n")
  cat("Number of units:",sum(s),"\n")
  cat("Beta coefficient(s):", r$coefficients,"\n")
  cat("Std. error:", r$std_error,"\n")
  cat("t-value:", r$t_value, "\n")
  cat("p_value:",r$p_value,"\n")  
  cat("cov_matrix:\n")
  print(r$cov_matrix)
    }  
  }
if(description) 
  cat("The regression estimator is:\n")
s1
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/rhg.r"
rhg<-function(X,selection)
{
if(missing(selection)) stop("'selection' is missing")
if(!is.data.frame(X)) X=as.data.frame(X)
if(ncol(X)<2) stop("the input is incorrect")
if(any(is.na(X))) stop("the input should not contain missing values")
if(is.null(names(X))) stop("specify the names of columns in X")
if(!all(c("ID_unit","status") %in% names(X))) stop("ID_unit and/or status are not columns of X") 
m=match(selection,names(X),nomatch = 0)
if(sum(m)==0) stop("'selection' should be the name of one of the X's columns")
X1=cbind.data.frame(X$status,X[,m])
names(X1)=c("status",names(X)[m])
x=unique(X1[,2])
if(ncol(X1)>=3)
for(i in 3:ncol(X1))
  x=list(x,unique(X1[,i]))
x=expand.grid(x)
ng=1
prob=rhgroup=numeric(nrow(X1))
for (i in 1:nrow(x)) {
            expr=rep(FALSE, nrow(X1))
            for(j in 1:nrow(X1)) {
                                  expr[j] = all(X1[j,2:ncol(X1)] == x[i, ])
                                  if(expr[j]) rhgroup[j]=ng
                                 }
            if(any(expr)) ng=ng+1 
                     }
gr=unique(rhgroup)
if(is.data.frame(X1))
X1=cbind.data.frame(X1,rhgroup)
else X1=cbind(X1,rhgroup)
for(i in 1:length(gr))
  {l=nrow(X1[X1[,ncol(X1)]==gr[i],]) 
   lr=nrow(X1[X1[,ncol(X1)]==gr[i] & X1[,1]==1,])
   for(j in 1:length(prob)) 
             if(rhgroup[j]==gr[i] & X1[j,1]==1) prob[j]=lr/l
   }
result=cbind.data.frame(X$ID_unit,X1,prob)
names(result)=c("ID_unit",names(X1),"prob_resp")
res = NULL
mm = match(names(X), names(result), nomatch = 0)
if(0 %in% mm)
{index = (1:ncol(X))[mm == 0]
res = cbind.data.frame(X[X$ID_unit==result$ID_unit, index], result)
names(res)[1:length(index)] = names(X)[index]
}
else res=result
res
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/rhg_strata.r"
rhg_strata<-function(X,selection)
{
if(is.matrix(X)) X=as.data.frame(X)
m=match(selection,names(X),nomatch=0)
if(sum(m)==0) stop("the 'selection' should be the name of one the X columns")
if(!("Stratum" %in% names(X))) stop("the column 'Stratum' is missing")
result=NULL
u=unique(X$Stratum)
for(i in 1:length(u))
  {si=X[X$Stratum==u[i],]
   x=cbind.data.frame(si$ID_unit,si$status,si[,m])
   names(x)=c("ID_unit","status",names(X)[m])
   result=rbind.data.frame(result,rhg(x,selection))
   }
res = NULL
mm = match(names(X), names(result), nomatch = 0)
index = (1:ncol(X))[mm == 0]
if (length(index) > 0) {
res = cbind.data.frame(X[X$ID_unit==result$ID_unit, index], result)
names(res)[1:length(index)] = names(X)[index]
                        }
res
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/rmodel.r"
rmodel<-function(formula,weights,X)
{
cl <- match.call() 
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "weights"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- model.response(mf, "numeric")
w <- as.vector(model.weights(mf))
x <- model.matrix(mt, mf, contrasts)
prob<-glm(y~x,family="binomial",weights=w)$fitted.values
result<-cbind.data.frame(X,prob)
names(result)<-c(names(X),"prob_resp")
result
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/samplecube.R"
"samplecube" <-
function(X,pik,order=1,comment=TRUE,method=1) 
{
EPS = 1e-11
N=length(pik)
if(is.array(X)!=TRUE) X=array(X,c(N,length(X)/N)) 
if(method==1)
{
    if (length(pik[pik > EPS & pik < (1 - EPS)]) > 0) 
        pikstar = fastflightcube(X, pik, order, comment)
    else {
        if (comment == TRUE) 
            cat("\nNO FLIGHT PHASE")
        pikstar = pik
    }
    if (length(pikstar[pikstar > EPS & pikstar < (1 - EPS)]) > 
        0) 
        pikfin = landingcube(X, pikstar, pik, comment)
    else {
        if (comment == TRUE) 
            cat("\nNO LANDING PHASE")
        pikfin = pikstar
    }
}
else
{   p=length(X)/length(pik)
    pikstar=pik
    for(i in 0:(p-1))
      {
      if (length(pikstar[pikstar > EPS & pikstar < (1 - EPS)]) > 0) 
      pikstar = fastflightcube(X[,1:(p-i)]/pik*pikstar, pikstar, order, comment)
      }
      pikfin = pikstar
      for(i in 1:N) if(runif(1)<pikfin[i]) pikfin[i]=1
}
if (comment == TRUE) {
        A <- X[pik > EPS, ]/pik[pik > EPS]
        TOT = t(A) %*% pik[pik > EPS]
        EST = t(A) %*% pikfin[pik > EPS]
        DEV = 100 * (EST - TOT)/TOT
        cat("\n\nQUALITY OF BALANCING\n")
        if(is.null(colnames(X)))  Vn = as.character(1:length(TOT)) else Vn=colnames(X);
        for(i in 1:length(TOT)) if(Vn[i]=="") Vn[i]=as.character(i)
	d = data.frame(TOTALS = c(TOT), 
        HorvitzThompson_estimators = c(EST), Relative_deviation = c(DEV))
        rownames(d)<-Vn
        print(d)
    }
    round(pikfin)
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/srswor.R"
"srswor" <-
function(n,N)
{s<-rep(0,times=N);s[sample(N,n)]<-1;s}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/srswor1.R"
"srswor1" <-
function(n,N)
{j=0
 s=numeric(N)
 for(k in 1:N) if(runif(1)<(n-j)/(N-k+1)) {j=j+1;s[k]=1;}
 s
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/srswr.R"
"srswr" <-
function(n,N) as.vector(rmultinom(1,n,rep(n/N,times=N)))

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/strata.r"
strata<-function(data, stratanames=NULL, size, method=c("srswor","srswr","poisson","systematic"),pik,description=FALSE)
{
if(missing(method)) {warning("the method is not specified; by default, the method is srswor")
                     method="srswor"
                    }
if(!(method %in% c("srswor","srswr","poisson","systematic"))) 
  stop("the name of the method is wrong")
if(method %in% c("poisson","systematic") & missing(pik)) stop("the vector of probabilities is missing")
if(missing(stratanames)|is.null(stratanames))
   {
   if(method=="srswor")
	result=data.frame((1:nrow(data))[srswor(size,nrow(data))==1],rep(size/nrow(data),size))
   if(method=="srswr")
     {
      s=srswr(size,nrow(data))
      st=s[s!=0]
      l=length(st)
      result=data.frame((1:nrow(data))[s!=0])
      if(size<=nrow(data))    
          result=cbind.data.frame(result,st,prob=rep(size/nrow(data),l))
      else 
          {prob=rep(size/nrow(data),l)/sum(rep(size/nrow(data),l))
           result=cbind.data.frame(result,st,prob)
          }  
  
      colnames(result)=c("ID_unit","Replicates","Prob")
     }
   if(method=="poisson")
    {
     pikk=inclusionprobabilities(pik,size)
     s=(UPpoisson(pikk)==1)
     if(length(s)>0)
     result=data.frame((1:nrow(data))[s],pikk[s])
     if(description) 
      cat("\nPopulation total and number of selected units:",nrow(data),sum(s),"\n") 
     }
    if(method=="systematic")
    {
     pikk=inclusionprobabilities(pik,size)
     s=(UPsystematic(pikk)==1)
     result=data.frame((1:nrow(data))[s],pikk[s])
     }
   if(method!="srswr")
          colnames(result)=c("ID_unit","Prob")
   if(description & method!="poisson") cat("\nPopulation total and number of selected units:",nrow(data),sum(size),"\n")
   }
else
{
data=data.frame(data)
index=1:nrow(data)
m=match(stratanames,colnames(data))
if(any(is.na(m))) stop("the names of the strata are wrong")
data2=cbind.data.frame(data[,m])
x1=data.frame(unique(data[,m]))
colnames(x1)=stratanames
result=NULL
for(i in 1:nrow(x1))
	{
expr=rep(FALSE,nrow(data2))
for(j in 1:nrow(data2)) 
    expr[j]=all(data2[j,]==x1[i,])
y=index[expr]
if(description & method!="poisson") 
 {cat("Stratum" ,i,"\n")
  cat("\nPopulation total and number of selected units:",length(y),size[i],"\n")
  }
if(method!="srswr" & length(y)<size[i]) 
	{stop("not enough obervations for the stratum ", i, "\n") 
	 st=c(st,NULL)
       }
else
	{if(method=="srswor")
		{st=y[srswor(size[i],length(y))==1]
             r=x1[i,1]
             if(ncol(x1)>=2) 
             for(j in 2:ncol(x1))
                 r=cbind.data.frame(r,rep(x1[i,j],size[i]))
             r=cbind.data.frame(r,st,rep(size[i]/length(y),size[i]))
             }
       if(method=="systematic")
		{pikk=inclusionprobabilities(pik[y],size[i])
             s=(UPsystematic(pikk)==1)
             st=y[s]
             r=x1[i,1]
             if(ncol(x1)>=2)   
             for(j in 2:ncol(x1))
                 r=cbind.data.frame(r,rep(x1[i,j],size[i]))
             r=cbind.data.frame(r,st,pikk[s])
            }
 
	 if(method=="srswr") 
	   {s=srswr(size[i],length(y))
          st=cbind.data.frame(y[s!=0],s[s!=0])
          l=length(s[s!=0])
          r=x1[i,1]
          if(ncol(x1)>=2)  
          for(j in 2:ncol(x1))
             r=cbind.data.frame(r,rep(x1[i,j],l))   
          if(size[i]<=length(y))    
          r=cbind.data.frame(r,st,prob=rep(size[i]/length(y),l))
          else 
          {prob=rep(size[i]/length(y),l)/sum(rep(size[i]/length(y),l))
          r=cbind.data.frame(r,st,prob)}  
          }
       if(method=="poisson") 
		{pikk=inclusionprobabilities(pik[y],size[i])
             s=(UPpoisson(pikk)==1)
             if(any(s))
             { 
             st=y[s]
             r=x1[i,1]
             if(ncol(x1)>=2) 
 		 for(j in 2:ncol(x1))
                 r=cbind.data.frame(r,rep(x1[i,j],length(st)))
             r=cbind.data.frame(r,st,pikk[s])
            if(description) 
		 {cat("Stratum" ,i,"\n")
		  cat("\nPopulation total and number of selected units:",length(y),length(st),"\n")
		  }
               }
             else 
		{if(description) 
		 {cat("Stratum" ,i,"\n")
		  cat("\nPopulation total and number of selected units:",length(y),0,"\n")
              }
             r=NULL 
     	       }
            }
  
	}
if(!is.null(r))
	{
r=cbind(r,i)
result=rbind.data.frame(result,r)
        } 

 }  
if(method=="srswr")
colnames(result)=c(stratanames,"ID_unit","Replicates","Prob","Stratum")
else
colnames(result)=c(stratanames,"ID_unit","Prob","Stratum")
if(description)  {cat("Number of strata ",nrow(x1),"\n")
                  if(method=="poisson")
   			 cat("Total number of selected units", nrow(result),"\n")
                  else       
                  cat("Total number of selected units", sum(size),"\n")
                  }
}
result
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/varest.r"
varest<-function(Ys,Xs=NULL,pik,w=NULL)
{
 if (any(is.na(pik))) 
        stop("there are missing values in pik")
 if (any(is.na(Ys))) 
        stop("there are missing values in y")
 if (length(Ys) != length(pik)) 
        stop("y and pik have different sizes")
if(!is.null(Xs))
 {if(is.data.frame(Xs)) Xs=as.matrix(Xs)
  if(is.vector(Xs) & (length(Ys)!= length(Xs)))  
        stop("x and y have different sizes")
  if(is.matrix(Xs) & (length(Ys) != nrow(Xs)))  
        stop("x and y have different sizes")
}
 a=(1-pik)/sum(1-pik)
 if(is.null(Xs))
	{A=sum(a*Ys/pik)
	 var=sum((1-pik)*(Ys/pik-A)^2)/(1-sum(a^2))
	}
 else 
   {B=t(Xs*w) 
   beta=ginv(B%*%Xs)%*%B%*%Ys
   e=Ys-Xs%*%beta
   A=sum(a*e/pik)
   var=sum((1-pik)*(e/pik-A)^2)/(1-sum(a^2))
}
var
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/vartaylor_ratio.r"
vartaylor_ratio=function(Ys,Xs,pikls)
{
if (any(is.na(pikls))) 
        stop("there are missing values in pikls")
if(nrow(pikls)!=ncol(pikls))
        stop("pikls is not a square matrix")
if (any(is.na(Ys))) 
        stop("there are missing values in y")
if (any(is.na(Xs))) 
        stop("there are missing values in x")
if (length(Ys) != nrow(pikls) | length(Xs) != nrow(pikls) | length(Xs) != length(Ys)) 
      stop("y, x and pikls have different sizes")
pik=diag(pikls)
n=length(pik)
xhat=sum(Xs/pik)
yhat=sum(Ys/pik)
r=yhat/xhat
z=(Ys-r*Xs)/xhat
delta=matrix(0,nrow=nrow(pikls),ncol=ncol(pikls))
for(i in 1:(n-1))
 {for(j in (i+1):n)
      delta[i,j]=delta[j,i]=(1-pik[i]*pik[j]/pikls[i,j])*z[i]*z[j]/(pik[i]*pik[j])
  delta[i,i]=(1-pik[i])*z[i]^2/pik[i]^2}
delta[n,n]=(1-pik[n])*z[n]^2/pik[n]^2
list(ratio=r, estvar=sum(delta))
}



#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/writesample.R"
"writesample" <-
function(n,N) 
{
if(n==N) samples=rep(1,times=N) 
else{ 
x=numeric(N) 
row=1 
for(i in (n+1):N) row=row*i 
k=1 
for(i in 1:(N-n)) k=k*i 
row=row/k 
samples=matrix(0,row,N) 
k=1 
sol=0 
x[1]=-1 
while(k<=N && k>0) 
{ 
while(x[k]<1) 
{ 
x[k]=x[k]+1 
s=0 
for(i in 1:N) 
s=s+x[i] 
if(s==n) { 
sol=sol+1 
samples[sol,]=x 
} 
else 
if(k<N){k=k+1 
x[k]=-1 
} 
} 
k=k-1 
} 
} 
samples
}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/RtmpXIUJTn/R.INSTALL10d63af1/sampling/R/zzz.R"
.First.lib <-function (lib, pkg) {
library.dynam("sampling", pkg, lib)
}
