src/MinPakke.jl

@@ -44,11 +44,6 @@ export TRLooCV
 export PlotTRLooCV
 export TRLooCVNum
 export TRGCVNum
-export TRSegCV
-export TRVirCV
-export TRBidiagDecomp
-
-export simulateSpectrum
 
 
 include("preprocessing.jl")
@@ -56,6 +51,5 @@ include("convenience.jl")
 include("conveniencePlots.jl")
 include("variousRegressionFunctions.jl")
 include("TR.jl")
-include("simulateSpectroscopicData.jl")
 
 end

src/TR.jl

@@ -1,6 +1,13 @@
 
-using Optimization        # For numerical minimization of PRESS statistic
+
-using OptimizationOptimJL # For numerical minimization of PRESS statistic
+
+
+
+
+using Optimization
+using OptimizationOptimJL
+
+
 
 struct TRSVD
     U::Matrix{Float64}
@@ -15,26 +22,12 @@ struct TRSVD
 end
 
 
-struct TRBidiag
-    W::Matrix{Float64}
-    B::Bidiagonal{Float64, Vector{Float64}}
-    T::Matrix{Float64}
-    y::Vector{Float64}
-    mX::Matrix{Float64}
-    my::Float64
-    regType::String
-    regParam1::Float64
-    regMat::Matrix{Float64}
-    n::Int64
-    p::Int64
-end
-
 
 """
 ### TO DO: ADD FRACTIONAL DERIVATIVE REGULARIZATION ###
 
-    regularizationMatrix(X; regType="L2", regParam1=0, regParam2=1e-14)
+    regularizationMatrix(X; regType="legendre", regParam1=0, regParam2=1e-14)
-    regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
+    regularizationMatrix(p::Int64; regType="legendre", regParam1=0, regParam2=1e-14)
 
 Calculates and returns square regularization matrix.
 
@@ -62,56 +55,22 @@ end
 
 function regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
     
-if regType == "bc" # Discrete derivative with boundary conditions
+if regType == "bc"
     regMat  = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
-elseif regType == "legendre" # Fill in polynomials in bottom row(s) to get square matrix
+elseif regType == "legendre"
     regMat  = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
     P, _    = plegendre(regParam1-1, p);   
     regMat[end-regParam1+1:end,:] = sqrt(regParam2) * P;
 elseif regType == "L2"
     regMat = I(p);
-elseif regType == "std" # Standardization
+elseif regType == "std"
     regMat = regParam2;
-elseif regType == "GL" # Grünwald-Letnikov fractional derivative regulariztion
-    # regParam1 is alpha (order of fractional derivative)
-    C = ones(p)*1.0;
-    for k in 2:p
-        C[k] = (1-(regParam1+1)/(k-1)) * C[k-1];
-    end
-
-    regMat = zeros(p,p);
-
-    for i in 1:p
-        regMat[i:end, i] = C[1:end-i+1];
-    end
 end
 
 return regMat
 end
 
 
-"""
-    function TRBidiagDecomp(X, y, A=(minimum(size(X))-1), regType="L2", regParam1=0, regParam2=1e-14)
-
-Calculates regularization matrix (using function "RegularizationMatrix"),
-and centres and transforms data matrix according to "X / regMat".
-Output is an object of type "TRBidiag" and is used as input to other TR functions.
-"""
-function TRBidiagDecomp(X, y, A=(minimum(size(X))-1), regType="L2", regParam1=0, regParam2=1e-14)
-
-n, p       = size(X);
-mX         = mean(X, dims=1);    
-X          = X .- mX;
-my         = mean(y);
-y          = vec(y .- my);
-regMat     = regularizationMatrix(X; regType, regParam1, regParam2);
-X          = X / regMat;
-_, W, B, T = bidiag2(X, y, A);
-TRObj      = TRBidiag(W, B, T, y, mX, my, regType, regParam1, regMat, n, p);
-    
-return TRObj
-end
-
 
 """
     function TRSVDDecomp(X, regType="L2", regParam1=0, regParam2=1e-14)
@@ -133,112 +92,6 @@ TRObj   = TRSVD(U, s, V, mX, regType, regParam1, regMat, n, p);
 return TRObj
 end
 
-function TRSegmentOrth(X, segments)
-
-n = size(X,1);
-n_segments = maximum(segments);
-
-U = zeros(n,n);
-
-for seg in 1:n_segments
-    inds = vec(seg .== segments)
-    U[inds, inds], _, _ = svd(X[inds,:], full=false);
-end
-
-return U
-end
-
-"""
-    function TRVirCV(X, y, segments, lambdas, regType="L2", regParam1=0, regParam2=1e-14)
-
-Segmented virtual cross-validation (VirCV) for TR models.
-Outputs: b, press, lambda_min, lambda_min_ind, GCV
-b are (virtual) press-minimal regression coefficients.
-"""
-function TRVirCV(X, y, segments, lambdas, regType="L2", regParam1=0, regParam2=1e-14)
-
-U_segments = TRSegmentOrth(X, segments);
-bs         = vec(sum(U_segments, dims=1).^2);
-
-n, p    = size(X);
-mX      = mean(X, dims=1);    
-X       = X .- mX;
-my      = mean(y);
-y       = vec(y .- my);
-X       = U_segments' * X;
-y       = U_segments' * y;
-regMat  = regularizationMatrix(p; regType, regParam1, regParam2);
-X       = X / regMat;
-U, s, V = svd(X, full=false);
-
-denom  = broadcast(.+, broadcast(./, lambdas, s'), s')';
-H      = broadcast(.+, U.^2 * broadcast(./, s, denom), bs./n);
-resid  = broadcast(.-, y, U * broadcast(./, s .* (U'*y), denom));
-rescv  = broadcast(./, resid, broadcast(.-, 1, H));
-press  = vec(sum(rescv.^2, dims=1));
-#rmsecv = sqrt.(1/n .* press);
-GCV    = vec(broadcast(./, sum(resid.^2, dims=1), mean(broadcast(.-, 1, H), dims=1).^2));
-
-lambda_min, lambda_min_ind = findmin(press);
-lambda_min_ind             = lambda_min_ind[1];
-
-denom2  = broadcast(.+, lambda_min ./ s', s')';
-b       = V * broadcast(./, (U' * y), denom2);
-b       = regMat \ b;
-b       = [my .- mX*b; b];
-
-return b, press, lambda_min, lambda_min_ind, GCV
-end
-
-
-
-"""
-    function TRSegCV(X, y, lambdas, cv, regType="L2", regParam1=0, regParam2=1e-14)
-
-Segmented cross-validation based on the Sherman-Morrison-Woodbury updating formula.
-Inputs:
-    - X : Data matrix
-    - y : Response vector
-    - lambdas : Vector of regularization parameter values
-    - cv : Vector of length n indicating segment membership for each sample
-    - regType, regParam1, regParam2 : Inputs to regularizationMatrix function
-
-Outputs: rmsecv, b, lambda_min, lambda_min_ind.
-b are regression coefficients corresponding to the lambda value minimising the CV-error.
-"""
-function TRSegCV(X, y, lambdas, cv, regType="L2", regParam1=0, regParam2=1e-14)
-
-TR        = TRSVDDecomp(X, regType, regParam1, regParam2);
-n_seg     = maximum(cv);
-n_lambdas = length(lambdas);
-my        = mean(y);
-y         = y .- my;
-
-denom  = broadcast(.+, broadcast(./, lambdas, TR.s'), TR.s')';
-resid  = broadcast(.-, y, TR.U * broadcast(./, TR.s .* (TR.U'*y), denom));
-rescv  = zeros(TR.n, n_lambdas);
-sdenom = sqrt.(broadcast(./, TR.s, denom))';
-
-for seg in 1:n_seg
-        
-    Useg = TR.U[vec(cv .== seg),:];
-    Id   = 1.0 * I(size(Useg,1)) .- 1/TR.n;
-        
-    for k in 1:n_lambdas
-        Uk                        = Useg .* sdenom[k,:]';
-        rescv[vec(cv .== seg), k] = (Id - Uk * Uk') \ resid[vec(cv .== seg), k];
-    end
-end
-
-press  = sum(rescv.^2, dims=1)';
-rmsecv = sqrt.(1/TR.n .* press);
-    
-lambda_min, lambda_min_ind = findmin(rmsecv)
-lambda_min_ind             = lambda_min_ind[1]
-b                          = TRRegCoeffs(TR, y, lambda_min, my)
-    
-return b, rmsecv, lambda_min, lambda_min_ind
-end
 
 
 """

src/simulateSpectroscopicData.jl

@@ -1,43 +0,0 @@
-Lk(waves, c, gamma) = @. gamma / (pi*(waves-c)^2+gamma^2);
-Gk(waves, c, gamma) = @. 1/(sqrt(2*pi)*gamma) * exp(-(waves-c)^2 / (2*gamma^2));
-
-function pseudoVoigtPeak(waves, c, gamma, eta, alpha)
-    
-fk = zeros(length(waves));
-
-for i in 1:length(c)
-        fk += alpha[i] .* pseudoVoigtPeak(waves, c[i], gamma[i], eta[i]);
-end
-    
-return fk
-end
-    
-function pseudoVoigtPeak(waves, c, gamma::Float64, eta::Float64)
-
-fk = eta * Lk(waves, c, gamma) + (1-eta) * Gk(waves, c, gamma);
-
-return fk
-end
-
-function simulateBaseline(waves, a)
-
-vm = zeros(length(waves), length(a));
-
-for i in 0:(length(a)-1)
-    vm[:, i+1] = waves.^i ./ norm(waves.^i);
-end
-
-b = vm * a;
-    
-return b
-end
-
-function simulateSpectrum(waves, c, gamma, eta, alpha, a, sigma=0.0)
-    
-pure_spec = pseudoVoigtPeak(waves, c, gamma, eta, alpha);
-b         = simulateBaseline(waves, a);
-noise     = randn(length(waves)) .* sigma;
-spec      = pure_spec + b + noise;
-return spec, pure_spec, b, noise
-end
-