


source("FstepArchetypesMod.R")
source("FstepArchetypoids.R")
source("FFunctions_to_calculate_real_archetypes_with_swap2.R")

library(archetypes)
library(fda)

temps <- matrix(scan("monthtemp.dat",0), 35, 12, byrow=T)


#  ----------------   names of weather stations  -------------------

meteonames <-   c("St. John_s",    "Charlottetown", "Halifax" ,
                  "Sydney",        "Yarmouth",      "Fredericton",
                  "Arvida",        "Montreal",      "Quebec City",
                  "Schefferville", "Sherbrooke",    "Kapuskasing",
                  "London",        "Ottawa",        "Thunder Bay",
                  "Toronto",       "Churchill",     "The Pas",
                  "Winnipeg",      "Prince Albert", "Regina",
                  "Beaverlodge",   "Calgary",       "Edmonton",
                  "Kamloops",      "Prince George", "Prince Rupert",
                  "Vancouver",     "Victoria",      "Dawson",
                  "Whitehorse",    "Frobisher Bay", "Inuvik",
                  "Resolute",      "Yellowknife")

#  ----------------   labels for months  ----------------------

monthletter <- c("J", "F", "M", "A", "M", "J",
                 "J", "A", "S", "O", "N", "D")

months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
            "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

#  ------------  indices for four meteorological regions  ----------

atlindex <- c(1,2,3,4,5,6,7,8,9,10,11,13,14,16)
pacindex <- c(25,26,27,28,29)
conindex <- c(12,15,17,18,19,20,21,22,23,24,30,31,35)
artindex <- c(32,33,34)



#  ---------------  argument vectors  --------------------

monthtime <- seq(0.5, 11.5, 1)    #  mid points of months


#  -------------------  set up the fourier basis  ----------------------

#    Note that an odd-number of basis functions are used, so that the
#      number of basis functions will be 13.  The functions can allow
#      for more basis functions than points being fit.

monthbasisfd <- create.fourier.basis( c(0,12), 12)

#  --------------  make temperature fd object  ------------

tempmat <- t(temps)
dimnames(tempmat) <- list(months, meteonames)


tempfd <- Data2fd(tempmat, monthtime, monthbasisfd,fdnames=c("Months", "Station", "Deg C"))





X=t(tempfd$coefs) #coefficeinets

PM=eval.penalty(monthbasisfd, rng=c(0,12))

K=4
set.seed(2015)

norep=20
lass10 <- FstepArchetypesMod(data=X,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes
 
an4 <-  FstepArchetypoids(K,TRUE,X,lass10,FALSE,K-1,PM)

aw4 <-  FstepArchetypoids(K,FALSE,X,lass10,FALSE,K-1,PM) #which_alpha archetypoids

ini_arch=c()
 for (j in 1:K){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
abeta <- FrealArchetypes(X,huge = 200, K, ini_arch,PM) #which beta archetypoid


# 8 34 30 29


##convert into functions
a4=ai$archetypes #son coeficientes en base de archetypes; archetypes

yaf=fd(t(a4),monthbasisfd)
plot(yaf)


#plot together
 
fdmat <- eval.fd(seq(0,12,12/500),yaf)
  matplot(seq(0,12,12/500), fdmat, type = "l",lty=1,xlab="Month",ylab="Temperature (deg. C)",cex.lab=1.5,cex.axis=1.5) #archetypes
  
fdmat2 <- eval.fd(seq(0,12,12/500),tempfd[c(8,34,30,29),])

lines(seq(0,12,12/500), fdmat2[,1], type = "l",lty=2,xlab="",ylab="",col=1) #archetypoids
lines(seq(0,12,12/500), fdmat2[,2], type = "l",lty=2,xlab="",ylab="",col=2)
lines(seq(0,12,12/500), fdmat2[,3], type = "l",lty=2,xlab="",ylab="",col=3)
lines(seq(0,12,12/500), fdmat2[,4], type = "l",lty=2,xlab="",ylab="",col=4)  


