using Jchemo, JchemoData
using JLD2, DataFrames, CairoMakie, GLMakie
using FreqTables

Data importation

path_jdat = dirname(dirname(pathof(JchemoData)))
db = joinpath(path_jdat, "data/challenge2018.jld2") 
@load db dat
@names dat

(:X, :Y)

X = dat.X
@head X

... (4075, 680)

Y = dat.Y
@head Y

... (4075, 4)

wlst = names(X) 
wl = parse.(Int, wlst)

680-element Vector{Int64}:
 1120
 1122
 1124
 1126
 1128
 1130
 1132
 1134
 1136
 1138
    ⋮
 2462
 2464
 2466
 2468
 2470
 2472
 2474
 2476
 2478

plotsp(X, wl; nsamp = 500, xlabel = "Wavelength (nm)").f

Preprocessing by SNV and derivation

model1 = snv()
model2 = savgol(npoint = 21, deriv = 2, degree = 3)
model = pip(model1, model2)
fit!(model, X)
Xp = transf(model, X)
@head Xp

3×680 Matrix{Float64}:
 -0.00393533  -0.00441755  -0.00477681  …  0.000514478  0.000481081
 -0.00121436  -0.0013095   -0.00135921     0.000996321  0.000930884
 -0.00355712  -0.00399785  -0.00434838     0.00044084   0.000421751
... (4075, 680)

plotsp(Xp, wl; nsamp = 500, xlabel = "Wavelength (nm)").f

Types of materials

typ = Y.typ
freqtable(string.(typ, " - ", Y.label))

10-element Named Vector{Int64}
Dim1                      │ 
──────────────────────────┼────
ANF - animal feed         │ 391
CLZ - rapeseed(ung)       │ 420
CNG - corn gluten         │ 395
EHH - grass silage        │ 422
FFS - full fat soya       │ 432
FRG - wheat (ung)         │ 411
MPW - milk powder & whey  │ 410
PEE - maize wp            │ 407
SFG - sun flower seed(gr) │ 281
TTS - soya meal           │ 506

test = Y.test  # training/test (0/1) observations
tab(test) 
freqtable(typ, test)

10×2 Named Matrix{Int64}
Dim1 ╲ Dim2 │   0    1
────────────┼─────────
ANF         │ 351   40
CLZ         │ 378   42
CNG         │ 356   39
EHH         │ 380   42
FFS         │ 397   35
FRG         │ 371   40
MPW         │ 372   38
PEE         │ 367   40
SFG         │ 272    9
TTS         │ 457   49

Split Train/Test from Tot

s = Bool.(test) # same as: s = Y.test .== 1
Xtrain = rmrow(Xp, s)
typtrain = rmrow(typ, s)
Xtest = Xp[s, :]
typtest = typ[s]
ntot = nro(X)
ntrain = nro(Xtrain)
ntest = nro(Xtest)
(ntot = ntot, ntrain, ntest)

(ntot = 4075, ntrain = 3701, ntest = 374)

Pca on Train

nlv = 3
model = pcasvd(; nlv)
fit!(model, Xtrain)
@head T = model.fitm.T

3×3 Matrix{Float64}:
 -0.00444931   0.0171604    0.00192674
 -0.0119228   -0.00286706   0.00500189
 -0.00391203   0.0180218   -0.000975336
... (3701, 3)

CairoMakie.activate!()  
#GLMakie.activate!()    # for interactive axe-rotation

i = 1
plotxyz(T[:, i], T[:, i + 1], T[:, i + 2]; size = (600, 500), color = (:red, .3), markersize = 10, 
    xlabel = string("PC", i), ylabel = string("PC", i + 1), zlabel = string("PC", i + 2), 
    title = "Pca score space").f

lev = mlev(typtrain)
nlev = length(lev)
colm = cgrad(:tab10, nlev; categorical = true, alpha = .5)
i = 1
plotxyz(T[:, i], T[:, i + 1], T[:, i + 2], typtrain; size = (700, 500), color = colm, markersize = 10, 
    xlabel = string("PC", i), ylabel = string("PC", i + 1), zlabel = string("PC", i + 2), 
    title = "Pca score space").f

Projection of the test observations

@head Ttest = transf(model, Xtest)

3×3 Matrix{Float64}:
  0.0126087    0.0102876    0.00597325
 -0.00690901  -0.00507066   0.00472288
 -0.00284267   0.0236669   -7.55497e-5
... (374, 3)

f, ax = plotxyz(T[:, i], T[:, i + 1], T[:, i + 2], typtrain; size = (700, 500), color = colm, markersize = 10, 
    xlabel = string("PC", i), ylabel = string("PC", i + 1), zlabel = string("PC", i + 2), 
    title = "Pca score space")
scatter!(ax, Ttest[:, i], Ttest[:, i + 1], Ttest[:, i + 2]; markersize = 6, color = :black)
f

---