

source("LstepArchetypesMod.R")
source("LstepArchetypoids.R")
source("LFunctions_to_calculate_real_archetypes_with_swap2.R")

library(archetypes)

d=read.table('responses.txt',header=T)

##################
#ADA
X=d[,2:22]

nada=10
norep=20
set.seed(2018)
adi=stepLArchetypoids3(X, 1:nada, norep=norep)

screeplot(adi)

adi[[3]]
#[[1]]
#[1] 182 274   1

apply(X[c(182,274,1),],1,sum)/21
 #     182       274         1 
#0.5714286 0.5238095 0.0952381


#ternary plot


a3=X[adi[[3]][[1]],]
data=X
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')

 a=alphascce

colnames(a)=c('A1','A2','A3')

library(ggtern)
ggtern(data=data.frame(a),aes(x=A1,y=A2,z=A3),aes(x,y,z))  +
stat_density_tern(
geom='polygon',
aes(fill=..level..),
bins=45, base='identity',
color='grey')  +
scale_fill_gradient(low='yellow',high='red')+
 limit_tern(1.01,1.01,1.01)+
geom_point()


 wa=apply(a,1,which.max)
table(wa)

####################


#AA



norep=20
set.seed(2018)
aa3=LstepArchetypesMod(X, 3, nrep=norep)

ai=bestModel(aa3)

ai$archetypes
round(ai$archetypes, digits = 2) 
round(ai$archetypes>.5, digits = 2) 

 ba=ai$archetypes>.5
 apply(ba,1,sum)/21
#[1] 0.0952381 0.4761905 0.6190476



####################
#cluster PAM

library(cluster)


 da=daisy(X,type = list(symm =1:21))

pa=pam(da,k=3)

table(pa$clustering)

 


n=690
p=princomp(da)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:(n),1],pr[2]*v[1:(n),2])

plot(pv,xlab='Dim 1',ylab='Dim 2',type='n')
points(pv[pa$clustering==1,], col=1,pch=1)
points(pv[pa$clustering==2,], col=2,pch=2)
points(pv[pa$clustering==3,], col=4,pch=3)



############################

#homals

library(homals)

rh3 <- homals(X,,ndim=3,level=rep("nominal",21))
rh3
plot3dstatic(rh3, plot.type = "loadplot")
plot3dstatic(rh3, plot.type = "objplot")

rh <- homals(X,level=rep("nominal",21))
rh
plot(rh)
plot(rh,plot.type="objplot")
plot(rh,plot.type="objplot",main='')

##############################
#LCA cluster


library(poLCA)


# Note for poLCA coding must start at 1.
d1=d+1
#d2=apply(d1[,2:22],2,cut,2,labels=c(1,2))

#attach(d2)

#Create the formula object for poLCA command 
f = cbind(p1  , p2  , p3,   p4 ,  p5,   p6a , p6b , p7  , p8 ,  p9 ,  p10 , p11 , p12 , p13a, p13b, p14a, p14b, p15,  p16,  p17a, p17b)~1

set.seed(1234)
#Run LCA for number of classes, nclass, from 1 to 5
lca1 = poLCA(f, data.frame(d1), nclass = 1)
lca2 = poLCA(f, data.frame(d1), nclass = 2)
lca3 = poLCA(f, data.frame(d1), nclass = 3, maxiter = 3500)
lca4 = poLCA(f, data.frame(d1), nclass = 4, maxiter = 21000)
lca5 = poLCA(f, data.frame(d1), nclass = 5, maxiter = 21000)


#Create an elbow plot for the number of classes and the BIC values
plot(1:5, c(lca1$bic, lca2$bic, lca3$bic, lca4$bic,lca5$bic), type = "b",
     pch = 19, xlab = "Number of Clusters", ylab = "BIC", main = "", xaxt = "n")
axis(1, 1:5, 1:5)

#To look at the outcome for the 3-class fit
lca3

#plot(rh$objscores, col=lca3$predclass)


plot(rh$objscores,xlab='Dimension 1',ylab='Dimension 2',type='n')
points(rh$objscores[lca3$predclass==1,], col=1,pch=1)
points(rh$objscores[lca3$predclass==2,], col=2,pch=2)
points(rh$objscores[lca3$predclass==3,], col=4,pch=3)

prof=matrix(0,nrow=21,ncol=3)
for (i in 1: 21)
{prof[i,]=as.vector(lca3$probs[[i]][,2])
}

round(t(prof),2)
apply(round(t(prof)>.5,2),1,sum)/21

table(lca3$predclass)

######################
#k-means

set.seed(1234)
k3=kmeans(d[,2:22],3,nstart=10)

apply(k3$centers>.5,1,sum)/21
k3$centers

round(k3$centers,2)

table(k3$cluster)

####################
#PAA in MATLAB


%data

 d=dlmread('responses.txt','',1,0);

X=d(:,2:22);

%code based on compareAA_cluster Seth & Eugster, PAMI 2016 

nSamFeat=X;

rng default, flagSave = false;
epsilon = 10^-12; maxTrial = 1; 
computeLikelihood = @(x,p)(sum(x .* log(p + epsilon) + (1-x) .* log(1 - p + epsilon)));



% Number of prototypes
        nLat = 3;
        
        % Construct input for AA
        clear nFeatSam_
        for countFeat = 1:size(nSamFeat,2)
            nFeatSam_{countFeat}(1,:) = nSamFeat(:,countFeat);
            nFeatSam_{countFeat}(2,:) = 1 - nSamFeat(:,countFeat);
        end
        
        clear matSamLat matLatSam objAA
        options = generate_options();
        for countTrial = 1:maxTrial
            [matSamLat{countTrial}, matLatSam{countTrial}, objAA{countTrial}] = paa_nominal_EM(nFeatSam_, nLat, options);
        end
        % Best AA solution indexed by indAA
        [~, indAA] = max(cellfun(@(x)(x(end)),objAA));
        % Archetypes
        matLatFeatAA = (nSamFeat' * matSamLat{indAA})';
        
 ######################
#Summary

s=apply(X,1,sum)*100/21  

ecdf(s)(57.1)

ecdf(s)(52.38)

ecdf(s)(33.3)

 ecdf(s)(47.6)
 ecdf(s)(61.9)

ecdf(s)(38.1)
ecdf(s)(42.9)


quantile(s, probs = seq(0, 1, 0.01))

#############################


library(e1071)
dh=hamming.distance(as.matrix(X))

 dh[182,274]
dh[586,162]

hamming.distance(ba[2,],ba[3,])