anchor.diff:                      else 
                       matrix (caPrev |> List.rev)
               loop 0 []
  -
  +        ///calculate null space of Ca (based on SVD)
           let Z =
               if Ca.NumRows * Ca.NumCols = 0 then 
                   Matrix.identity n            
               else
  -                let k = 
  +                let nullspace = 
                       let (eps,U,V) =
                           Algebra.LinearAlgebra.SVD Ca
                       let rank = eps |> Seq.filter (fun x -> x >= 1e-5) |> Seq.length //Threshold if a singular value is considered as 0. //mat.NumRows - eps.Length
                       let count = V.NumRows - rank
                       Matrix.getRows V (rank) (count)
                       |> Matrix.transpose
  -                k
  +                nullspace
          
           let I = Matrix.identity n
   
           let G_lambda (lambd:float)  = 
               let dInverse = Algebra.LinearAlgebra.Inverse D
  -            2. * ( H.Transpose * dInverse * H + lambd / (float n)  * (W.Transpose * W))
  +            2. * ( H.Transpose * dInverse * H + lambd / (float n)  * W.Transpose * W)
       
           let A_tild (lambd:float)    = 
               let zproInverse = Algebra.LinearAlgebra.Inverse (Z.Transpose * (G_lambda lambd) * Z)
  -            2. * (lambd / (float n)) * ( Z * zproInverse) * Z.Transpose * (W.Transpose * W)
  -
  -        let V_tild (lambd:float) =     
  -            let aMat = a |> Matrix.ofVector
  -            let yMat = y |> Matrix.ofVector
  -            let no = 
  -                (W * (aMat - yMat))
  -                |> fun m -> Matrix.getCol m  0
  -                |> Vector.norm
  -                |> fun x -> x*x
  +            2. * (lambd / (float n)) * Z * zproInverse * Z.Transpose * W.Transpose * W
  +
  +        let V_tild rho (lambd:float) =     
  +            let no =
  +                W * (a - y)
  +                |> Vector.fold (fun acc x -> acc + (pown x 2)) 0.   
     
  -            let tr = I - (A_tild lambd) |> getDiag |> Vector.sum
  +            let tr = 
  +                I - rho * (A_tild lambd)
  +                |> Matrix.trace
   
               (float n) * no / (tr * tr)
   
           let var_ny (lambd:float) =
               let a = ((W * (a-y)) |>Vector.norm ) ** 2.
  -            let b = I - (A_tild lambd) |> getDiag |> Vector.sum
  +            let b = I - (A_tild lambd) |> Matrix.getDiag |> Vector.sum
               a/b
   
  -        let gcv = V_tild lambda
  +        //standard gcv
  +        let gcv  = V_tild 1.0 lambda
  +        //modified gcv
  +        let mGcv = V_tild 1.3 lambda
  +        //robust gcv
  +        let rGcv() =
  +            let A = A_tild lambda
  +            let gamma = 0.2
  +            let mu = ((A * A) |> Matrix.trace) / (float n)
  +            (gamma + (1. - gamma) * mu) * gcv
  +
   
           let variance = var_ny lambda
  -        (gcv,variance)
  -
  -    let private spline operator (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con=
  +        //(gcv,variance)
  +        (mGcv,gcv)
  +
  +    //record type for easy handling inside the shape classification process
  +    type InnerResult = {
  +        TraceA  : Vector<float>
  +        GCV     : float
  +        Var     : float
  +        Error   : float
  +        CTemp   : float [,]
  +        Lambda  : float}
  +
  +    let createInnerResult a g c e v l = {
  +        TraceA  = a
  +        GCV     = g
  +        Var     = v
  +        Error   = e
  +        CTemp   = c
  +        Lambda  = l}
  +
  +    ///calculates a spline with the shape defined by the operator(increasing or decreasing) and the condition(0,1,2,or 3 extrema)
  +    let private spline operator (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con =
   
           /// Sets the diagonal to value inplace
           let setDiagonalInplace value (m:Matrix<_>) =
  
positive.code:                             Matrix.setRow Ctemp (4*i-1) (- (temp - (B.Row (i) |> vector)))
                              if i = j - 1 then
                                  Matrix.setRow Ctemp (4*i-4) (Vector.zero Ctemp.NumCols)
                                  Matrix.setRow Ctemp (4*i-3) -(B.Row (i-1) |> vector) //0??
                                  Matrix.setRow Ctemp (4*i-2) (Vector.zero Ctemp.NumCols)
                                  Matrix.setRow Ctemp (4*i-1) (Vector.zero Ctemp.NumCols)
                          for i = j to n-1 do
  
negative.code:                 diagVecNew |> Vector.map (fun x -> x / trace * length)
              for i = 0 to xVal.Length - 1 do Wtemp.[i,i] <- finalDiag.[i]
              Wtemp
          | CV -> 
              let Wtemp = Matrix.create xVal.Length xVal.Length 0.
              let cvOfVec =
                  if yVal.Length%numRep = 0 then
  