anchor.diff:              let gcvfinal = gcvmin.[eminIndex]
               let lambdafinal = lambdafin.[eminIndex]
               let cfinal = cfin.[eminIndex]
  -            createHermiteResultSimple afinal cfinal gcvfinal efinal lambdafinal ctemp nan
  -
  -    let splineIncreasing (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con= 
  +            let splineFunction = initEvalAt x afinal cfinal
  +            createTempClassResult afinal cfinal gcvfinal efinal lambdafinal ctemp infinity splineFunction
  +
  +
  +
  +    ///calculates a constrained initially increasing spline with given y-,and y-Values, a weighting matrix and the number of allowed extrema
  +    let splineIncreasing (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con = 
           spline (~-) x y W con
  -
  -    let splineDecreasing (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con= 
  +    
  +    ///calculates a constrained initially decreasing spline with given y-,and y-Values, a weighting matrix and the number of allowed extrema
  +    let splineDecreasing (x:Vector<float>) (y:Vector<float>) (W:Matrix<float>) con = 
           spline (~+) x y W con
  -        
  -    
  -    let getBestFitOfWeighting arr wMat xVal=
  -        let getinitialestimate = getInitialEstimateOfWXY wMat (vector arr) xVal 
  -        let initial = getinitialestimate |> fun (a,b,c,d) -> b
  -        
  -        let fst1 =(1,splineIncreasing xVal (vector arr) wMat 0 )  |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.04) x.Lambda x.Ctemp x.Variance)
  -        let fst2 =(-1,splineDecreasing xVal (vector arr) wMat 0)  |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.04) x.Lambda x.Ctemp x.Variance)
  -        let bestfst = [fst1;fst2] |> List.minBy (fun (cl,x) -> x.Error)                                                                                                    
  -                                                                                                                                                                            
  -        let snd1 =(2,splineIncreasing xVal (vector arr) wMat 1; ) |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.06) x.Lambda x.Ctemp x.Variance)
  -        let snd2 =( 2,splineDecreasing xVal (vector arr) wMat 1;) |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.06) x.Lambda x.Ctemp x.Variance)
  -        let bestsnd = [snd1;snd2] |> List.minBy (fun (cl,x) -> x.Error)                                                                                                    
  -                                                                                                                                                                             
  -        let trd1 =(3,splineIncreasing xVal (vector arr) wMat 2; ) |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.08) x.Lambda x.Ctemp x.Variance)
  -        let trd2 =(-3,splineDecreasing xVal (vector arr) wMat 2;) |> fun (cl,x) -> (cl, createHermiteResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.08) x.Lambda x.Ctemp x.Variance)
  -        let besttrd = [trd1;trd2] |> List.minBy (fun (cl,x) -> x.Error)
  -        
  +
  +    ///returns the best fit of the observation x,-and y-values and a given weightingMatrix
  +    let getBestFitOfWeighting x_values y_values wMat= 
  +        let getinitialestimate = getInitialEstimateOfWXY wMat (vector y_values) x_values
  +        let initialSelectionCriterion = getinitialestimate.AICc
  +        
  +        //monotone splines
  +        let fst1 =(In0,splineIncreasing x_values (vector y_values) wMat 0) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.04) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let fst2 =(De0,splineDecreasing x_values (vector y_values) wMat 0) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.04) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let bestfst = [fst1;fst2] |> List.minBy (fun (cl,x) -> x.AICc)                                                                
  +        //splines with one extremum
  +        let snd1 =(In1,splineIncreasing x_values (vector y_values) wMat 1) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.06) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let snd2 =(De1,splineDecreasing x_values (vector y_values) wMat 1) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.06) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let bestsnd = [snd1;snd2] |> List.minBy (fun (cl,x) -> x.AICc)                                                               
  +        //splines with two extrema                                                                                                                                             
  +        let trd1 =(In2,splineIncreasing x_values (vector y_values) wMat 2) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.08) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let trd2 =(De2,splineDecreasing x_values (vector y_values) wMat 2) |> fun (shapeClass,(constraintMatrices,result)) -> shapeClass,result //(cl, createTempClassResultSimple x.TraceA x.TraceC x.Error (x.GCV*1.08) x.Lambda x.Ctemp x.AICc x.SplineFunction)
  +        let besttrd = [trd1;trd2] |> List.minBy (fun (cl,x) -> x.AICc)
  +
  +        //selection of optimal shape possibility by model selection via AICc 
           [bestfst;bestsnd;besttrd]
  -        |> List.minBy (fun (cl,x) -> x.GCV)
  -        |> fun (cl,x) -> if x.GCV < initial then (cl,x) else 0,x
  -
  -    ///gets indices of minima and maxima
  +        |> List.minBy (fun (cl,result) -> result.AICc)
  +        |> fun (cl,result) -> if result.AICc < initialSelectionCriterion then (cl,result) else Complex,result        
  +
  +        
  +
  +    ///gets indices of minima and maxima by brute force search
       let private investigateTrace pct (arr:float[]) =
           let almostEqual (pct:float) f1 f2 =
               let (min,max) = f2 - Math.Abs(f2 * pct),f2 + Math.Abs(f2 * pct)
  
positive.code:                             Matrix.setRow Ctemp (4*i-1) ((temp - (B.Row (i) |> vector)))
                          Matrix.setRow Ctemp (4*(n-1)) (B.Row (n-1) |> vector)
                          list <- (Matrix.copy Ctemp)::list
                  //printfn "ctemp fertig"
                  list
                  |> Array.ofList 
  
              else failwith "Condition max 2"
  
  
          let getError (y:Vector<float>) (a:Vector<float>) (W:Matrix<float>) =
              let tmp =  W * (y - a)
              tmp |> Seq.averageBy (fun x -> x**2.) |> sqrt
  
      
          let n = x.Length
  
          //Define intervals    (4*(n-1)+1) n
          let h = Array.init (n-1) (fun i -> x.[i+1] - x.[i] )
          
          //generation of matrices D and H (Wood)
          let H = Array2D.zeroCreate (n-2) n
          let D = Array2D.zeroCreate (n-2) (n-2)
  
  
negative.code:                     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 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
                  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)
      
          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 rho (lambd:float) =     
              let no =
                  W * (a - y)
                  |> Vector.fold (fun acc x -> acc + (pown x 2)) 0.   
    
              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) |> Matrix.getDiag |> Vector.sum
              a/b
  
          //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)
          (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<_>) =
  