[chronojump] Implementing crossvalidate for splines



commit bd869d2a34e7afe5e3da7abdc0d9b6f232d32081
Author: Xavier de Blas <xaviblas gmail com>
Date:   Sun Apr 17 12:19:32 2016 +0200

    Implementing crossvalidate for splines

 encoder/call_capture.R     |    2 +
 encoder/call_graph.R       |    2 +
 encoder/graph.R            |   35 ++++++++++++++++------
 encoder/graphSmoothingEC.R |   28 ++++++++++++++++--
 encoder/util.R             |   66 +++++++++++++++++++++++++++++++++++++++++--
 src/encoder.cs             |    7 +++-
 src/encoderCapture.cs      |    8 +++++
 src/encoderRProc.cs        |    7 +++-
 src/gui/encoder.cs         |    3 ++
 src/utilEncoder.cs         |    2 +-
 10 files changed, 138 insertions(+), 22 deletions(-)
---
diff --git a/encoder/call_capture.R b/encoder/call_capture.R
index 5889b36..a1c088e 100644
--- a/encoder/call_capture.R
+++ b/encoder/call_capture.R
@@ -27,6 +27,8 @@ source(scriptCaptureR)
 DEBUG <- FALSE
 DebugFileName <- paste(options[5], "/chronojump-debug.txt", sep="")
 
+CROSSVALIDATESMOOTH <- FALSE
+
 
 f <- file("stdin")
 open(f)
diff --git a/encoder/call_graph.R b/encoder/call_graph.R
index 1690b59..caebb63 100644
--- a/encoder/call_graph.R
+++ b/encoder/call_graph.R
@@ -31,6 +31,8 @@ Translated = unlist(strsplit(options[29], "\\;"))
 DEBUG <- FALSE
 DebugFileName <- paste(options[5], "/chronojump-debug.txt", sep="")
 
+CROSSVALIDATESMOOTH <- FALSE
+
 source(paste(options[4], "/util.R", sep=""))
 source(paste(options[4], "/graphSmoothingEC.R", sep=""))
 
diff --git a/encoder/graph.R b/encoder/graph.R
index 5be533f..5708159 100644
--- a/encoder/graph.R
+++ b/encoder/graph.R
@@ -2138,6 +2138,7 @@ doProcess <- function(options)
        unicodeWorks <<- checkUnicodeWorks()
 
        DEBUG <<- op$Debug
+       CROSSVALIDATESMOOTH <<- op$CrossValidate
        
        print(c("1 Title=",op$Title))
        #unicoded titles arrive here like this "\\", convert to "\", as this is difficult, do like this:
@@ -2392,9 +2393,14 @@ doProcess <- function(options)
                
                #print(curves, stderr())
        
-               #find SmoothingsEC
-               SmoothingsEC = findSmoothingsEC(singleFile, displacement, curves, op$Eccon, op$SmoothingOneC,
-                                               NULL, NULL, NULL, NULL) #this row is only needed for 
singleFile (signal)
+               #find SmoothingsEC. TODO: fix this
+               if(CROSSVALIDATESMOOTH) {
+                       for(i in 1:n)
+                               SmoothingsEC[i] = 0
+               }
+               else
+                       SmoothingsEC = findSmoothingsEC(singleFile, displacement, curves, op$Eccon, 
op$SmoothingOneC,
+                                                       NULL, NULL, NULL, NULL) #this row is only needed for 
singleFile (signal)
                print(c("SmoothingsEC:",SmoothingsEC))
        } else {        #singleFile == True. reads a signal file
                displacement=scan(file=op$File,sep=",")
@@ -2468,10 +2474,16 @@ doProcess <- function(options)
                #print(curves)
                
                #find SmoothingsEC
-               SmoothingsEC = findSmoothingsEC(
-                                               singleFile, displacement, curves, op$Eccon, op$SmoothingOneC,
-                                               op$EncoderConfigurationName, op$diameter, op$inertiaMomentum, 
op$gearedDown 
-                                               ) #second row is needed for singleFile (signal)
+               if(CROSSVALIDATESMOOTH) {
+                       for(i in 1:n)
+                               SmoothingsEC[i] = 0
+               }
+               else
+                       SmoothingsEC = findSmoothingsEC(
+                                                       singleFile, displacement, curves, op$Eccon, 
op$SmoothingOneC,
+                                                       op$EncoderConfigurationName, op$diameter, 
op$inertiaMomentum, op$gearedDown 
+                                                       ) #second row is needed for singleFile (signal)
+
                print(c("SmoothingsEC:",SmoothingsEC))
                
                if(curvesPlot) {
@@ -2661,14 +2673,17 @@ doProcess <- function(options)
                                debugParameters(listN(xUpperValue, xLowerValue), "paint all smoothing 2")
 
                                #3.a) find max power (y) for some smoothings(x) (closer)
-                               x <- seq(from = xUpperValue, to = xLowerValue, length.out = 5)
+                               x <- seq(from = xUpperValue, to = xLowerValue, length.out = 8)
                                y <- smoothAllSetYPoints(x, displacementAllSet, 
                                                         op$EncoderConfigurationName, op$MassBody, 
op$MassExtra, op$ExercisePercentBodyWeight, 
                                                         op$gearedDown, op$anglePush, op$angleWeight, 
op$diameter, op$inertiaMomentum)
 
                                #3.b) create a model with x,y to find optimal x (in closer values)
-                               smodel <- smooth.spline(y,x)
-                               smoothingAll <- predict(smodel, maxPowerAtAnyRep)$y
+                               #but only if we have 4+ unique values for smooth.spline. If not, previous 
smoothingAll will be used
+                               if(length(unique(x)) >= 4 && length(unique(y)) >= 4) {
+                                       smodel <- smooth.spline(y,x)
+                                       smoothingAll <- predict(smodel, maxPowerAtAnyRep)$y
+                               }
                        }
 
                        debugParameters(listN(x, y, maxPowerAtAnyRep, smoothingAll), "paint all smoothing 3")
diff --git a/encoder/graphSmoothingEC.R b/encoder/graphSmoothingEC.R
index 0ba4d2d..9596994 100644
--- a/encoder/graphSmoothingEC.R
+++ b/encoder/graphSmoothingEC.R
@@ -157,13 +157,24 @@ findSmoothingsEC <- function(singleFile, displacement, curves, eccon, smoothingO
                                #2 get max power concentric (y) at eccentric-concentric phase with current 
smoothing of an interval of possible smoothings (x)
                                
                                x <- seq(from = as.numeric(smoothingOneC), 
-                                                   to = as.numeric(smoothingOneC)/2, 
-                                                   length.out = 5)
+                                                   to = as.numeric(smoothingOneC)/4, 
+                                                   length.out = 8)
                                y <- findSmoothingsECYPoints(eccentric.concentric, conStart, conEnd, x, 
maxPowerConAtCon, 
                                                myEncoderConfigurationName, myDiameter, 100, 
myInertiaMomentum, myGearedDown)
                                #write(paste("x, y", x, y), stderr())
 
 
+                               write("smooth.spline x (a)", stderr())  
+                               write(x, stderr())
+                               write("smooth.spline y (a)", stderr())  
+                               write(y, stderr())
+                       
+                               #if less than 4 unique x or y cannot smooth spline. Just use smoothingOneC as 
default value
+                               if(length(unique(x)) < 4 || length(unique(y)) < 4) {
+                                       smoothings[i] = smoothingOneC
+                                       next
+                               }
+
                                smodel <- smooth.spline(y,x)
                                smoothingOneEC <- predict(smodel, maxPowerConAtCon)$y
                                write(paste("smoothingOneEC", smoothingOneEC), stderr())
@@ -201,11 +212,22 @@ findSmoothingsEC <- function(singleFile, displacement, curves, eccon, smoothingO
                                x <- seq(
                                         from = xUpperValue, 
                                         to = xLowerValue,
-                                        length.out = 5)
+                                        length.out = 8)
                                y <- findSmoothingsECYPoints(eccentric.concentric, conStart, conEnd, x, 
maxPowerConAtCon, 
                                                myEncoderConfigurationName, myDiameter, 100, 
myInertiaMomentum, myGearedDown)
                                
+                       
+                               write("smooth.spline x (b)", stderr())
+                               write(x, stderr())
+                               write("smooth.spline y (b)", stderr())  
+                               write(y, stderr())
                                
+                               #if less than 4 unique x or y cannot smooth spline. Just use recent calculed 
smoothingOneEC as default value
+                               if(length(unique(x)) < 4 || length(unique(y)) < 4) {
+                                       smoothings[i] = smoothingOneEC
+                                       next
+                               }
+
                                smodel <- smooth.spline(y,x)
                                smoothingOneEC <- predict(smodel, maxPowerConAtCon)$y
                                        
diff --git a/encoder/util.R b/encoder/util.R
index c2a49ed..ac20346 100644
--- a/encoder/util.R
+++ b/encoder/util.R
@@ -72,7 +72,8 @@ assignOptions <- function(options) {
                    DecimalSeparator    = options[25],
                    Title               = options[26],
                    OperatingSystem     = options[27],  #if this changes, change it also at start of this R 
file
-                   Debug               = options[30]
+                   Debug               = options[30],
+                   CrossValidate       = options[31]
                    #Unassigned here:
                    #   englishWords [28]
                    #   translatedWords [29]
@@ -258,7 +259,12 @@ getSpeed <- function(displacement, smoothing) {
        #use position because this does not make erronously change the initial and end of the curve
        #do spline with displacement is buggy because data is too similar -2, -1, 0, 1, 2, ...
        position <- cumsum(displacement)
-       positionSpline <- smooth.spline( 1:length(position), position, spar=smoothing)
+       
+       spar <- smoothing
+       if(CROSSVALIDATESMOOTH)
+               spar <- cvParCall(1:length(position), position)
+
+       positionSpline <- smooth.spline( 1:length(position), position, spar=spar)
        
        #get speed converting from position to displacement, but repeat first value because it's missing on 
the diff process
        #do not hijack spline like this because this works for speed, but then acceleration has big error
@@ -271,14 +277,19 @@ getSpeed <- function(displacement, smoothing) {
        return (speedSpline)
 }
 
-#same smoothing than getSped
+#same smoothing than getSpeed
 getPositionSmoothed <- function(displacement, smoothing) {
        #no change affected by encoderConfiguration
        
        #use position because this does not make erronously change the initial and end of the curve
        #do spline with displacement is buggy because data is too similar -2, -1, 0, 1, 2, ...
        position <- cumsum(displacement)
-       positionSpline <- smooth.spline( 1:length(position), position, spar=smoothing)
+
+       spar <- smoothing
+       if(CROSSVALIDATESMOOTH)
+               spar <- cvParCall(1:length(displacement), displacement)
+
+       positionSpline <- smooth.spline( 1:length(position), position, spar=spar)
        
        return (positionSpline$y)
 }
@@ -328,6 +339,7 @@ reduceCurveBySpeed <- function(eccon, startT, startH, displacement, smoothingOne
        speed <- getSpeed(displacement, smoothing)
        
        speed.ext=extrema(speed$y)
+       print(speed.ext)
 
        #in order to reduce curve by speed, we search the cross of speed (in 0m/s)
         #before and after the peak value, but in "ec" and "ce" there are two peak values:
@@ -1115,7 +1127,53 @@ getInertialDiametersPerMs <- function(displacement, diametersPerTick)
   return(diameter)
 }
 
+#----------- Begin spar with crossvalidation (Aleix Ruiz de Villa) -------------
+
+cvParCall <- function(x, y) {
+       parRange <- seq(0.1, 1, 0.05)
+       cvResults <- cvPar(x, y, parRange)
+
+       ## Agafem la spar que te un error més petit
+       #print( paste('The best spar is: ', parRange[ which.min(cvResults) ] ) )
+       write("CrossValidation:", stderr())
+       write(parRange[which.min(cvResults)], stderr())
+       return (parRange[which.min(cvResults)])
+}
+cvPar <- function(x, y, parRange, cvProp = 0.8) {
+
+       n <- length(x)
+       nTrain <- floor( cvProp * n )
+       nTest <- n - nTrain
+
+       cvParError <- vector(length=length(parRange))
+       for( parPos in 1:length(parRange) ){
+               selectedPar <- parRange[ parPos ]
+
+               ## Separem la mostra en train i test
+               indexTrain <- sample(1:n, nTrain)
+               indexTest <- (1:n)[-indexTrain]
+               xTrain <- x[ indexTrain ]
+               xTest <- x[ indexTest ]
+               yTrain <- y[ indexTrain ]
+               yTest <- y[ indexTest ]
+
+               ## Fem servir les dades del train per a construir els splines
+               splineModel <- smooth.spline(xTrain, yTrain, spar = selectedPar)
+
+               ## Fem la prediccio sobre la mostra test
+               predictedBySpline <- predict(splineModel, xTest)
+               # plot(xTest, predictedBySpline$y, col = 'blue')
+               # lines(xTest, yTest, col = 'green', type = 'p')
+
+               ## Calculem els errors
+               predictionError <- sqrt( mean( (predictedBySpline$y - yTest)^2 ) )
+               cvParError[ parPos ] <- predictionError
+       }
+
+       return(cvParError)
+}
 
+#----------- end spar with crossvalidation -------------
 
 #----------- Begin debug file output -------------
 debugParameters <- function (parameterList, currentFunction) 
diff --git a/src/encoder.cs b/src/encoder.cs
index 1fbf2a5..abaea7d 100644
--- a/src/encoder.cs
+++ b/src/encoder.cs
@@ -165,6 +165,7 @@ public class EncoderGraphROptions
        public string englishWords;
        public string translatedWords;
        public bool debug;
+       public bool crossValidate;
        
        public EncoderGraphROptions(
                        string inputData, string outputGraph, string outputData1, 
@@ -172,7 +173,7 @@ public class EncoderGraphROptions
                        EncoderParams ep,
                        string title, string operatingSystem,
                        string englishWords, string translatedWords,
-                       bool debug)
+                       bool debug, bool crossValidate)
        {
                this.inputData = inputData;
                this.outputGraph = outputGraph;
@@ -185,6 +186,7 @@ public class EncoderGraphROptions
                this.englishWords = englishWords;
                this.translatedWords = translatedWords;
                this.debug = debug;
+               this.crossValidate = crossValidate;
        }
 
        public override string ToString() {
@@ -199,7 +201,8 @@ public class EncoderGraphROptions
                        "#operatingsystem\n" +  operatingSystem + "\n" +
                        "#englishWords\n" +     englishWords + "\n" +
                        "#translatedWords\n" +  translatedWords + "\n" +
-                       "#debug\n" +            Util.BoolToRBool(debug) + "\n";
+                       "#debug\n" +            Util.BoolToRBool(debug) + "\n" +
+                       "#crossValidate\n" +    Util.BoolToRBool(crossValidate) + "\n";
        }
        
 
diff --git a/src/encoderCapture.cs b/src/encoderCapture.cs
index ae924ac..e376e8b 100644
--- a/src/encoderCapture.cs
+++ b/src/encoderCapture.cs
@@ -319,6 +319,14 @@ public abstract class EncoderCapture
                                        if(startFrame < 0)
                                                startFrame = 0;
 
+                                       LogB.Information("TTTT - i," + i.ToString() +
+                                                       "; directionChangeCount: " + 
+                                                       directionChangeCount.ToString() + 
+                                                       "; lastNonZero: " +
+                                                       lastNonZero.ToString() +
+                                                       "; final: " + 
+                                                       ((i - directionChangeCount + 
lastNonZero)/2).ToString());
+
                                        ecc = new EncoderCaptureCurve(
                                                        startFrame,
                                                        (i - directionChangeCount + lastNonZero)/2      
//endFrame
diff --git a/src/encoderRProc.cs b/src/encoderRProc.cs
index 40a2183..ca5d812 100644
--- a/src/encoderRProc.cs
+++ b/src/encoderRProc.cs
@@ -31,6 +31,7 @@ public abstract class EncoderRProc
        public enum Status { WAITING, RUNNING, DONE } 
        public Status status;
        public bool Debug = false;
+       public bool CrossValidate;
 
        protected string optionsFile;   
        protected EncoderStruct es;
@@ -303,7 +304,8 @@ public class EncoderRProcCapture : EncoderRProc
                                es, 
                                false,  //neuromuscularProfile
                                false,  //translate (graphs)
-                               Debug
+                               Debug,
+                               false   //crossValidate (unactive on capture at the moment)
                                ).ToString();
 
                TextWriter writer = File.CreateText(optionsFile);
@@ -538,7 +540,8 @@ public class EncoderRProcAnalyze : EncoderRProc
                                es, 
                                neuromuscularProfileDo,
                                translate,
-                               Debug
+                               Debug,
+                               CrossValidate
                                ).ToString();
 
                TextWriter writer = File.CreateText(optionsFile);
diff --git a/src/gui/encoder.cs b/src/gui/encoder.cs
index 9fde2de..c2ee144 100644
--- a/src/gui/encoder.cs
+++ b/src/gui/encoder.cs
@@ -126,6 +126,7 @@ public partial class ChronoJumpWindow
        [Widget] Gtk.CheckButton check_encoder_analyze_show_force;
        [Widget] Gtk.CheckButton check_encoder_analyze_show_power;
        
+       [Widget] Gtk.CheckButton checkbutton_crossvalidate;
        [Widget] Gtk.Button button_encoder_analyze;
        [Widget] Gtk.Box hbox_encoder_analyze_progress;
        [Widget] Gtk.Button button_encoder_analyze_cancel;
@@ -4608,6 +4609,8 @@ public partial class ChronoJumpWindow
 
                        encoder_pulsebar_analyze.Text = Catalog.GetString("Please, wait.");
                        encoderRProcAnalyze.status = EncoderRProc.Status.WAITING;
+       
+                       encoderRProcAnalyze.CrossValidate = checkbutton_crossvalidate.Active;
                
                        encoderThread = new Thread(new ThreadStart(encoderDoAnalyze));
                        GLib.Idle.Add (new GLib.IdleHandler (pulseGTKEncoderAnalyze));
diff --git a/src/utilEncoder.cs b/src/utilEncoder.cs
index bdb5589..b822460 100644
--- a/src/utilEncoder.cs
+++ b/src/utilEncoder.cs
@@ -351,7 +351,7 @@ public class UtilEncoder
                                title, operatingSystem,
                                Util.StringArrayToString(Constants.EncoderEnglishWords,";"),
                                Util.StringArrayToString(encoderTranslatedWordsOK,";"),
-                               debug
+                               debug, crossValidate
                                );
        }
 


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]