PartlyOrdered=function(X,y,dist,hdenom,b0,n,phi00num,maxit,tol)
{
#input
# X matrix of covariates
# y response vector
# dist matrix with Riemannian distances between manifold observations
# hdenom=denominator of the bandwith of the kernel
# b0 initial values for the parameter b of the linear part 
# n number of observations
# phi00num initial value for the phi0 function
# maxit maxim number of iterations
# tol tolerance
#
#return
# res a list with the estimations of parameter b and phi0 function
 
 
phi00=rep(phi00num,n)

  # ----------
  # kernel
  # ----------
 
  h=pi/hdenom
  theta=matrix(rep(1,n*n),ncol=n)
  theta[dist>0.0001]=(abs(sin(dist[dist>0.0001]))/dist[dist>0.0001])^(m*(k-1/2-m/2)-2)
  disth=dist/h
  kernelVar = apply(disth, c(1,2), dnorm )
  w=kernelVar/theta
  denom = apply(w, 1, sum )
  w=w/denom


  # ----------
  # REGRESIÓN ORDINAL
  # ----------

  le=levels(y)
  p1=as.numeric(y==le[1])
  p2=as.numeric(y==le[2])
  p3=as.numeric(y==le[3])
  yreal=c(p1[1],p2[1]+p1[1])

  for(j in 2: n){
  ytemp=c(p1[j],p2[j]+p1[j])
  yreal=c(yreal,ytemp)
   }

  bi=b0
  bf=b0+100
  i=0
  bn=b0
  phi0n=phi00
  phi=w%*%X


  Xbar=matrix(c(X[1,],X[1,]),byrow=T, ncol=2)
  phibar=matrix(c(phi[1,],phi[1,]), byrow=T,ncol=2)

  for(j in 2: n){
    xn=matrix(c(X[j,],X[j,]), byrow=T, ncol=2)
    phin=matrix(c(phi[j,],phi[j,]), byrow=T, ncol=2)
    Xbar=rbind(Xbar,xn)
    phibar=rbind(phibar,phin)
  }

  ###La matriz del kernel S2###
  S2=matrix(rep(0,4*n*n),ncol=2*n)
  cont=1
  for(j in 1:n ){
       cont2=1
       for(l in 1:n){
        S2[cont:(cont+1),cont2:(cont2+1)]=diag(rep(w[j,l],2))
        cont2=cont2+2
        }
      cont=cont+2
     }

  #######################################
  #Iterations
  #######################################

  i=0
  while(i<=maxit & (norm(as.matrix(bi-bf),'f')/norm(as.matrix(bi),'f'))>tol){
    i=i+1
    bi=bn
    gn=phi0n-phibar%*%bn
    Pn=1/(1+exp(-Xbar%*%bn-gn))
    p1n=Pn[seq(1,2*n-1,2)]
    p2n=Pn[seq(2,2*n,2)]-p1n
    p3n=1-(p1n+p2n)
   
    ### W weight matrix
    ### D matrix with the derivative of the link function in the diagonal

    W=matrix(rep(0,4*n*n),ncol=2*n)
    D=matrix(rep(0,4*n*n),ncol=2*n)

    cont=1
    for(j in 1:n ){
       temp=matrix((1/p2n[j])*c((1-p3n[j])/p1n[j], -1, -1, (1- p1n[j])/p3n[j]),ncol=2)
       W[cont:(cont+1),cont:(cont+1)]=temp
       temp=matrix(c( p1n[j]*(1-p1n[j]), 0, 0, p3n[j]*(1-p3n[j])),ncol=2)
       D[cont:(cont+1),cont:(cont+1)]=temp
       cont=cont+2
       }
    Z=Xbar%*%bn+gn+D%*%(yreal-Pn)
    phi0n=S2%*%Z
    bn=solve((t(Xbar-phibar)) %*% W %*% (Xbar-phibar)) %*% (t((Xbar-phibar))) %*% W %*% (Z-phi0n)
    bf=bn
    }
    res=list(bn=bn,phi0n=phi0n)
    return(res)
   }
