anchor.diff:      open FSharp.Stats
       open FSharp.Stats.Optimization
   
  +    let private getDiag (m:Matrix<'a>)=
  +        let range = min m.NumRows m.NumCols
  +        vector [for i = 0 to range - 1 do yield m.[i,i]]
  +
  +    let unzip4 arr  =
  +        Array.foldBack (fun (a,b,c,d) (accA,accB,accC,accD) -> a::accA, b::accB, c::accC, d::accD) arr ([],[],[],[])
  +        |> fun (ra,rb,rc,rd) -> Array.ofList ra,Array.ofList rb,Array.ofList rc,Array.ofList rd
  +
  +    ///works if MKL service is available
  +    let getGCV (C:Matrix<float>) (D:Matrix<float>) (H:Matrix<float>) (W:Matrix<float>) (n:int) (y:vector) (a:vector) (lambda:float) =
  +        let constr = C*(a |> Matrix.ofVector) 
  +        let tol = 0.001 
  +        let Ca =    
  +            let rec loop i caPrev =
  +                if i < C.NumRows then
  +                    if System.Math.Abs(constr.[i,0]) <= tol then 
  +                        loop (i+1) (((C.Row i) |> Seq.toArray)::caPrev)
  +                    else loop (i+1) caPrev
  +                else 
  +                    if caPrev |> List.isEmpty then
  +                        Matrix.create 0 0 0.
  +                    else 
  +                    matrix (caPrev |> List.rev)
  +            loop 0 []
  +
  +            //let mutable Caprev :float[][] = [||]
  +            //for i = 0 to C.NumRows - 1 do
  +            //    if System.Math.Abs(constr.[i,0]) <= tol then 
  +            //        Caprev <- Array.append Caprev [|(C.Row i) |> Seq.toArray|]
  +
  +            //if Caprev = [||] then
  +            //    Matrix.create 0 0 0.
  +            //else matrix Caprev
  +       
  +        let Z =
  +            if Ca.NumRows * Ca.NumCols = 0 then 
  +                Matrix.identity n            
  +            else
  +                let k = 
  +                    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
  +       
  +        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))
  +    
  +        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
  +  
  +            let tr = I - (A_tild lambd) |> getDiag |> Vector.sum
  +
  +            (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
  +            a/b
  +
  +        let gcv = V_tild lambda
  +
  +        let variance = var_ny lambda
  +        (gcv,variance)
   
       let private spline operator (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) lambd =
   
  
positive.code:         let Q = calcGlambda (Matrix.ofArray2D D) (Matrix.ofArray2D H) W n lambd 
          let c = calcclambda W n y lambd |> Vector.toArray
  
          let b = Array.create n (-infinity,0.) // Reconsider
  
          let a' = 
              let tmp = Matrix.create Q.NumRows Q.NumCols 1. |> setDiagonalInplace 0.5
  
negative.code:             // padding zero
              Vector.init (tmp.Length+2) (fun i -> if i>0 && i<=tmp.Length then tmp.[i-1] else 0.0) 
  
          //works if MKL service either is available 
          //let (gcv',var') = getGCV (A |> Matrix.ofArray2D) (Matrix.ofArray2D D) (Matrix.ofArray2D H) W y.Length y a' lambd
          (a',e',b',c')
  
  
  