<|start_filename|>src/BioFSharp.Mz/SearchEngineResult.fs<|end_filename|>
﻿namespace BioFSharp.Mz

open BioFSharp
open SearchDB


module SearchEngineResult = 
    
    type SearchEngine = 
        | AndromedaLike 
        | SEQUESTLike   
        | XTandemLike

    type SearchEngineResult<'a> = {
        SearchEngine          : SearchEngine
        SpectrumID            : string;     
        ModSequenceID         : int
        PepSequenceID         : int;
        GlobalMod             : int;
        IsTarget              : bool;
        ScanTime              : float; 
        StringSequence        : string;
        PrecursorCharge       : int;
        PrecursorMZ           : float
        MeasuredMass          : float;
        TheoMass              : float;            
        PeptideLength         : int;
        Score                 : 'a;
        NormDeltaBestToRest   : float;
        NormDeltaNext         : float;}
    
    ///
    let createSearchEngineResult searchEngine spectrumID modSequenceID pepSequenceID globalMod isTarget scanTime stringSequence precursorCharge precursorMZ measuredMass theoMass peptideLength xcorr normDeltaBestToRest normDeltaNext = 
        { SearchEngine=searchEngine;SpectrumID = spectrumID; ModSequenceID = modSequenceID; PepSequenceID = pepSequenceID;ScanTime = scanTime ; GlobalMod = globalMod; IsTarget = isTarget; StringSequence = stringSequence; 
          PrecursorCharge = precursorCharge; PrecursorMZ=precursorMZ;  MeasuredMass = measuredMass; TheoMass = theoMass; PeptideLength = peptideLength; Score = xcorr; NormDeltaBestToRest = normDeltaBestToRest;NormDeltaNext = normDeltaNext }

    /// Calculates sequest-like delta normalized by best score to the best score.
    ///  (Xcorr(top hit) - Xcorr(n)) / Xcorr(top hit). Thus, the deltaCn for the top hit is
    ///  (Xcorr(top hit) - Xcorr(top hit)) / Xcorr(top hit) = 0.
    /// if the best Score equals 0. this function returns returns 1 for every PSM
    let calcNormDeltaBestToRest (sourceList:SearchEngineResult<'a> list) =
        match sourceList with
        | h1::rest -> 
            if h1.Score <= 0. then sourceList |> List.map (fun sls -> {sls with NormDeltaBestToRest = 1.})
            else
            sourceList
            |> List.map
                (fun sls ->
                    let normDeltaBestToRest  = (h1.Score - sls.Score) / h1.Score
                    { sls with NormDeltaBestToRest = normDeltaBestToRest } )      
        | []       -> []

    // Iterates over the score ranked PSMs and computes the score difference between adjacent
    // PSMs normalized by the Score of the best ranked PSM.
    /// if the best Score equals 0., this function returns returns 0 for every PSM.
    let calcNormDeltaNext (sourceList:SearchEngineResult<'a> list) =        
        let rec loop normF acc l = 
            match l with 
            | hLast::[] ->
                {hLast with NormDeltaNext = 0.}::acc
                |> List.rev 
            | hi::hii -> 
                let normDeltaNext = (hi.Score - hii.[0].Score) / normF 
                loop normF ({hi with NormDeltaNext = normDeltaNext}::acc) hii 
        match sourceList with
        | h1::rest -> 
            let normFactor = h1.Score
            if normFactor <= 0. then sourceList |> List.map (fun sls -> {sls with NormDeltaNext = 0.})
            else
            loop normFactor [] sourceList                  
        | []       -> []


<|start_filename|>src/BioFSharp.Mz/ProteinInference.fs<|end_filename|>
﻿namespace BioFSharp.Mz


open FSharpAux
open BioFSharp.PeptideClassification
open BioFSharp.IO.GFF3
open FSharpAux.IO.SchemaReader.Attribute
open System.Collections.Generic
open FSharp.Stats
open FSharp.Stats.Fitting.NonLinearRegression
open FDRControl
open FDRControl.MAYU


module ProteinInference = 
    open BioFSharp.Elements

    type internal StructuralEqualityComparer<'a when 'a : equality >() = 
        interface IEqualityComparer<'a> with
            member x.Equals(a,b) = a = b
            member x.GetHashCode a = hash a 
    
    /// For a single peptide Sequence, contains information about all proteins it might originate from and its evidence class.
    type ProteinClassItem<'sequence> = 
        {
            GroupOfProteinIDs: string []
            PeptideSequence: 'sequence
            Class: PeptideEvidenceClass
        }

    /// For a group of proteins, contains information about all peptides that might be used for its quantification.
    type InferredProteinClassItem<'sequence> = 
        {
            GroupOfProteinIDs: string
            PeptideSequence: 'sequence []
            Class: PeptideEvidenceClass
        }    

    /// For a group of proteins, contains information about all peptides that might be used for its quantification and score calculated for it.
    type InferredProteinClassItemScored =
        {
            GroupOfProteinIDs  : string
            PeptideSequence    : string []
            Class              : PeptideEvidenceClass
            TargetScore        : float
            DecoyScore         : float
            Decoy              : bool
            DecoyHasBetterScore: bool
            FoundInDB          : bool
        }

    /// For a group of proteins, contains information about all peptides that might be used for its quantification and score / q-value calculated for it.
    type InferredProteinClassItemQValue =
        {
            InfProtClassItem : InferredProteinClassItemScored
            QValue           : float
        }

    /// For a group of proteins, contains information about all peptides that are put into the output file.
    type Result =
        {
            GroupOfProteinIDs: string
            PeptideSequence  : string
            Class            : PeptideEvidenceClass
            TargetScore      : float
            DecoyScore       : float
            QValue           : float
        }

    type PSMInput =
        {
            [<FieldAttribute("PepSequenceID")>]
            PepSequenceID   : int
            [<FieldAttribute("StringSequence")>]
            Seq             :string
            [<FieldAttribute("PercolatorScore")>]
            Score : float
        }

    let private createInferredProteinClassItem proteinIDs evidenceClass peptideSequences = 
        {
            GroupOfProteinIDs = proteinIDs
            PeptideSequence = peptideSequences
            Class = evidenceClass
        }   

    let createProteinClassItem proteinIDs evidenceClass peptideSequence : ProteinClassItem<'sequence> = 
        {
            GroupOfProteinIDs = proteinIDs
            PeptideSequence = peptideSequence
            Class = evidenceClass
        }

    let createInferredProteinClassItemScored proteinIDs evidenceClass peptideSequences targetScore decoyScore isDecoy decoyHasBetterScore foundInDB =
        {
            GroupOfProteinIDs   = proteinIDs
            PeptideSequence     = peptideSequences
            Class               = evidenceClass
            TargetScore         = targetScore
            DecoyScore          = decoyScore
            Decoy               = isDecoy
            DecoyHasBetterScore = decoyHasBetterScore
            FoundInDB           = foundInDB
        }

    let createInferredProteinClassItemQValue infProtClassItemScored qValue=
        {
            InfProtClassItem = infProtClassItemScored
            QValue           = qValue
        }

    let createInferredProteinClassItemOut proteinIDs evidenceClass peptideSequences targetScore decoyScore qValue =
        {
            GroupOfProteinIDs = proteinIDs
            PeptideSequence   = peptideSequences
            Class             = evidenceClass
            TargetScore       = targetScore
            DecoyScore        = decoyScore
            QValue            = qValue
        }




 
    /// Used to decide wether overlapping groups of proteins should be kept or merged
    type IntegrationStrictness = 
        /// Results in the minimal set of proteins which still can explain all measured peptides
        | Minimal
        /// All protein groups stay intact
        | Maximal

    /// Used to decide which peptides should be used for quantification of protein groups
    type PeptideUsageForQuantification =
        /// Use only the best matching peptides
        | Minimal
        /// Use all available peptides which point to a protein group which includes the given protein group
        | Maximal
        /// Use all available peptides which point to a part of the protein group
        | MaximalInverse


    let private increaseClass c =
        match c with
        | PeptideEvidenceClass.C1a -> PeptideEvidenceClass.C1b
        | PeptideEvidenceClass.C1b -> PeptideEvidenceClass.C2a
        | PeptideEvidenceClass.C2a -> PeptideEvidenceClass.C2b
        | PeptideEvidenceClass.C2b -> PeptideEvidenceClass.C3a
        | PeptideEvidenceClass.C3a -> PeptideEvidenceClass.C3b

    let private mapClassName c =
        match c with
        | PeptideEvidenceClass.C1a -> "Class1A"
        | PeptideEvidenceClass.C1b -> "Class1B"
        | PeptideEvidenceClass.C2a -> "Class2A"
        | PeptideEvidenceClass.C2b -> "Class2B"
        | PeptideEvidenceClass.C3a -> "Class3A"
        | PeptideEvidenceClass.C3b -> "Class3B"

    /// Appends a collection of proteins to a single string
    let proteinGroupToString (proteinGroup:string[]) =
        Array.reduce (fun x y ->  x + ";" + y) proteinGroup

    let removeModification pepSeq =
        String.filter (fun c -> System.Char.IsLower c |> not && c <> '[' && c <> ']') pepSeq

    /// Checks if GFF line describes gene
    let isGene (item: GFFLine<seq<char>>) =
        match item with
        | GFFEntryLine x -> x.Feature = "gene"
        | _ -> false

    /// Checks if GFF line describes rna
    let isRNA (item: GFFLine<seq<char>>) =
        match item with
        | GFFEntryLine x -> if x.Feature = "mRNA" then Some x else None
        | _ -> None

    /// Reads geographical information about protein from gff entry and builds the modelinfo of it
    /// This function takes an RNA gff3 entry and therefore will contain the splice variant id of the gene in its result.
    /// This splice variant id should be the same in the given FastA-file.
    let createProteinModelInfoFromEntry i locus (entry:GFFEntry) =
        let attributes = entry.Attributes
        /// Same as in FastA file
        let spliceVariantID =
            match Map.tryFind "Name" attributes with
            | Some res ->
                res.Head
            | None ->
                failwithf "could not find spliceVariantId for locus %s" locus
        let chromosomeID = entry.Seqid

        let direction =
            match entry.Strand with
            |'+' -> StrandDirection.Forward
            |'-' -> StrandDirection.Reverse

        createProteinModelInfo spliceVariantID chromosomeID direction locus i Seq.empty Seq.empty

    /// By reading GFF creates the protein models (relationships of proteins to each other) which basically means grouping the rnas over the gene loci
    /// TODO: Don't group over order but rather group over id
    let assignTranscriptsToGenes regexPattern (gffLines: seq<GFFLine<seq<char>>>)  =
        gffLines
        // transcripts are grouped by the gene they originate from
        |> Seq.groupWhen isGene
        |> Seq.map (fun group ->
            match Seq.head group with
            | GFFEntryLine x ->
                let locus = x.Attributes.["ID"].Head // =genename, this value is used to assign mRNAs of the same gene together
                group
                |> Seq.choose isRNA //the transcripts of the gene are chosen
                |> Seq.mapi (fun i element ->
                    // every transcript of gene gets its own number i and other info is collected from element and used for info of protein
                    let modelInfo = createProteinModelInfoFromEntry i locus element
                    let r = System.Text.RegularExpressions.Regex.Match(modelInfo.Id,regexPattern)
                    // the gff3 id has to be matched with the sequence in the fasta file. therefore the regexpattern is used
                    if r.Success then
                        r.Value,
                        modelInfo
                    else
                        failwithf "could not match gff3 entry id %s with regexpattern %s. Either gff3 file is corrupt or regexpattern is not correct" modelInfo.Id regexPattern
                )

            | _ -> Seq.empty
        )
        |> Seq.concat
        |> Map.ofSeq

    /// Creates a lookup data base to assign peptides to the proteins they are contained in
    let createPeptideProteinRelation (protModels:seq<ProteinModel<'id,'chromosomeId,'geneLocus,'sequence list> option>) =
        let ppRelation = BidirectionalDictionary<'sequence,ProteinModelInfo<'id,'chromosomeId,'geneLocus>>()
        protModels
        |> Seq.iter (fun prot ->
            // insert peptide-protein relationship
            // Todo: change type of proteinID in digest
            match prot with
            | Some proteinModel ->
                proteinModel.Sequence
                |> Seq.iter (fun pepSequence -> ppRelation.Add pepSequence proteinModel.ProteinModelInfo)
            | None -> ()
        )
        ppRelation

    /// Contains all different classes with their associated proteins
    type private ClassMap =
        Map<PeptideEvidenceClass,BidirectionalDictionary<string,string>>

    module private ClassMap =
        ///Creates a new empty ClassMap
        let empty() : ClassMap =
            Map.empty
            |> Map.add PeptideEvidenceClass.C1a (BidirectionalDictionary<string,string>())
            |> Map.add PeptideEvidenceClass.C1b (BidirectionalDictionary<string,string>())
            |> Map.add PeptideEvidenceClass.C2a (BidirectionalDictionary<string,string>())
            |> Map.add PeptideEvidenceClass.C2b (BidirectionalDictionary<string,string>())
            |> Map.add PeptideEvidenceClass.C3a (BidirectionalDictionary<string,string>())
            |> Map.add PeptideEvidenceClass.C3b (BidirectionalDictionary<string,string>())

        ///Searches the Memory Map from highest ranking to lowest ranking class for any of the given proteinIds. If an overlap exists, returns its evidence class.
        let searchProts (m:ClassMap) (prots:string[]) =
            match Map.tryFindKey (fun c (d:BidirectionalDictionary<string,string>) ->
                    Array.exists (fun p -> d.ContainsKey p) prots
                  ) m with
            | Some c -> c
            | None -> PeptideEvidenceClass.Unknown

        ///Adds the protein group to the given class in the memory map
        let addGroup (m:ClassMap) (c:PeptideEvidenceClass) prots =
            let name = proteinGroupToString prots
            Seq.iter (fun p -> m.[c].Add p name) prots
            m

    // Creates a Map of peptides with their highest found score
    let createPeptideScoreMap (psmInputs: PSMInput list list) =
        psmInputs
        |> List.concat
        |> List.groupBy (fun psm -> psm.Seq)
        |> List.map (fun (sequence, psmList) ->
            // This sequence may contain modifications.
            // Depending on the type of lookup this map is used for, the modifications have to be removed.
            sequence,
            psmList
            |> List.maxBy (fun psm -> psm.Score)
            |> fun psm -> psm.Score
        )
        |> Map.ofList

    // Assigns a score to each protein with reverse digested peptides based on the peptides obtained in psm.
    let createReverseProteinScores (reverseProteins: (string*string[])[]) (peptideScoreMap: Map<string,float>) =
        // Remove modifications from map since protein inference was also done using unmodified peptides
        let scoreMapWithoutMods =
            peptideScoreMap
            |> Map.toArray
            |> Array.map (fun (seq, score) -> removeModification seq, score)
            |> Map.ofArray
        reverseProteins
        |> Array.map (fun (protein, peptides) ->
        protein,
            (
                peptides
                // looks wether the peptides resulting from the reverse digest appear in the peptides from the psm
                |> Array.map (fun pep ->
                    scoreMapWithoutMods.TryFind pep
                    |> (fun x ->
                        match x with
                        | Some score -> score
                        | None -> 0.
                    )
                )
                |> Array.sum,
                peptides
            )
        )
        // only hits are relevant
        |> Array.filter (fun (protein, (score, peptides)) -> score <> 0.)
        |> Map.ofArray

    /// Sums up score of all peptide sequences
    let assignPeptideScores (peptideSequences : string []) (peptideScoreMap : Map<string,float>) =
        peptideSequences
        |> Array.map (fun sequence -> peptideScoreMap.Item sequence)
        |> Array.sum

    /// Looks if the given protein accession is present in a map of identified decoy proteins and assigns its score when found.
    let assignDecoyScoreToTargetScore (proteins: string) (decoyScores: Map<string,(float * string[])>) =
        let prots = proteins |> String.split ';'
        prots
        |> Array.map (fun protein ->
            decoyScores.TryFind protein
            |> fun protOpt ->
                match protOpt with
                // peptides which pointed to the decoy version of this protein are discarded here, they can be included if needed
                | Some (score,peptides) -> score
                | None -> 0.
        )
        |> Array.max

    ///In a given class finds the overlap between the given and already added proteins and integrates them
    let private findAndIntegrate (m:ClassMap) (c:PeptideEvidenceClass) (proteins:string []) =
        let merge (dict:BidirectionalDictionary<string,string>) =
            let oldGroups =
                Array.choose (dict.TryGetByKey >> (Option.map Array.ofSeq)) proteins
                |> Array.concat
                |> Array.distinct
                |> Array.choose (fun pg -> Option.map (fun seqs -> seqs,pg) (dict.TryGetByValue pg))
            let newSet = Set.ofSeq proteins         
            oldGroups
            |> Array.iter (fun (x,y) -> printfn "%s" ((Seq.toArray >> proteinGroupToString) x) )
            oldGroups
            |> Array.iter (fun (seqs,group) ->
                let oldSet = Set.ofSeq seqs
                if Set.isProperSubset newSet oldSet then 
                    let newGroup = proteinGroupToString proteins
                    dict.RemoveValue group
                    Seq.iter (fun k -> dict.Add k newGroup) proteins
                elif Set.isProperSubset oldSet newSet then 
                    ()
                else
                    let protIntersect = Set.intersect newSet oldSet
                    let proteins = Set.toArray protIntersect
                    let newGroup = proteinGroupToString proteins
                    dict.RemoveValue group
                    Seq.iter (fun k -> dict.Add k newGroup) proteins
            ) 
            dict
        match ClassMap.searchProts m proteins with
        | PeptideEvidenceClass.Unknown ->
            printfn "Search did not return anything"
            ClassMap.addGroup m c proteins
        | foundC ->
            if foundC < c then
                m
            else
                Map.map (fun ca dict ->
                    if ca = c then merge dict
                    else dict) m
    
    let private createClassMap (integrationStrictness:IntegrationStrictness) (proteinClassItems: ProteinClassItem<'sequence> list ) =
        let rec loop (m:ClassMap) proteinClassItems remaining lookingForClass =
            match proteinClassItems with
            | (item:ProteinClassItem<'sequence>) :: l ->
                    if item.Class = lookingForClass then
                        let m' = 
                            if integrationStrictness = IntegrationStrictness.Maximal then 
                                ClassMap.addGroup m item.Class item.GroupOfProteinIDs
                            else  
                                (findAndIntegrate m item.Class item.GroupOfProteinIDs)
                        loop m' l remaining lookingForClass
                    else
                        loop m l (item::remaining) lookingForClass
            | [] ->
                if lookingForClass = PeptideEvidenceClass.C3b then
                    m
                else
                    loop m remaining [] (increaseClass lookingForClass)
        loop (ClassMap.empty()) proteinClassItems [] PeptideEvidenceClass.C1a

    /// Used to map the resulting protein groups to the sequences which are used for their quantification
    let createProteinToPepSequencesMap (proteinClassItems: ProteinClassItem<'sequence> list) =
        proteinClassItems
        |> List.collect (fun ci -> 
            let prots = ci.GroupOfProteinIDs
            List.init prots.Length (fun i -> prots.[i],ci)
           )
        |> List.groupBy fst
        |> List.map (fun (p,l) -> p, List.map snd l)
        |> Map.ofList
        //|> Map.ofList

    let inferSequences (integrationStrictness:IntegrationStrictness) (peptideUsageForQuantification:PeptideUsageForQuantification) (proteinClassItems: ProteinClassItem<'sequence> list) =
        ///Adds a value v to the set of values already associated with key k
        let classes : ClassMap = createClassMap integrationStrictness proteinClassItems
        /// Used to map the resulting protein groups to the sequences which are used for their quantification
        let seqDict = createProteinToPepSequencesMap proteinClassItems
        printfn "finish setup"
        classes
        |> Seq.collect (fun kv -> 
            let c = kv.Key
            let d = kv.Value
            d.GetArrayOfValues 
            |> Seq.map (fun group -> 
                let prots = Option.get (d.TryGetByValue group)
                let protSet = Set.ofSeq prots
                match peptideUsageForQuantification with
                | PeptideUsageForQuantification.Maximal ->                     
                    let potentialCandidates = 
                        prots
                        |> Seq.collect (fun prot -> 
                            seqDict.[prot]
                            |> Seq.choose (fun protClassItem -> 
                                let set = Set.ofArray protClassItem.GroupOfProteinIDs
                                if Set.isSubset protSet set then
                                    Some protClassItem
                                else 
                                    None
                                )
                            )                    
                    potentialCandidates
                    |> Seq.map (fun pci -> pci.PeptideSequence)
                    |> Seq.toArray
                    |> Array.distinct
                    |> createInferredProteinClassItem group c 
                | PeptideUsageForQuantification.Minimal -> 
                    let potentialCandidates = 
                        prots
                        |> Seq.collect (fun prot -> 
                            seqDict.[prot]
                            |> Seq.choose (fun protClassItem -> 
                                let set = Set.ofArray protClassItem.GroupOfProteinIDs
                                if Set.isSubset protSet set then
                                    Some protClassItem
                                else 
                                    None
                                )
                            )   
                    let withLength = 
                        potentialCandidates
                        |> Seq.map (fun pci -> pci.GroupOfProteinIDs.Length,pci)
                    let min = Seq.minBy fst withLength |> fst
                    Seq.filter (fst >> ((=) min)) withLength
                    |> Seq.map (snd >> (fun pci -> pci.PeptideSequence))
                    |> Seq.toArray
                    |> Array.distinct
                    |> createInferredProteinClassItem group c
                | PeptideUsageForQuantification.MaximalInverse ->    
                    let potentialCandidates = 
                        prots
                        |> Seq.collect (fun prot -> 
                            seqDict.[prot]
                            |> Seq.choose (fun protClassItem -> 
                                let set = Set.ofArray protClassItem.GroupOfProteinIDs
                                if Set.isSubset set protSet then
                                    Some protClassItem
                                else 
                                    None
                                )
                            )                    
                    potentialCandidates
                    |> Seq.map (fun pci -> pci.PeptideSequence)
                    |> Seq.toArray
                    |> Array.distinct
                    |> createInferredProteinClassItem group c                     
                )
            )

    module MAYU =

        /// Returns proteins sorted into bins according to their size.
        /// proteins are the proteins which were found in either the reverse or forward proteininference, proteinsFromDB are the proteins with peptide sequence
        /// on which the inference was performed.
        let binProteinsLength (proteins: InferredProteinClassItemScored []) (proteinsFromDB: (string*string)[]) binCount =
            // need to bin all proteins in the db, not only those with hit?
            let proteinsNotMatchedDB =
                let proteinsMatched =
                    proteins
                    |> Array.collect (fun protein ->
                        protein.GroupOfProteinIDs
                        |> String.split ';'
                    )
                    |> Set.ofArray
                // Creates an entry for every protein that is present in the search db and wasn't inferred
                proteinsFromDB
                |> Array.choose (fun (proteinName, peptideSequence) ->
                    if Set.contains proteinName proteinsMatched then
                        None
                    else
                        Some ((createInferredProteinClassItemScored proteinName BioFSharp.PeptideClassification.PeptideEvidenceClass.Unknown [|peptideSequence|] (-1.) (-1.) false false false),
                             float peptideSequence.Length)
                )

            // Adds the length of the peptide sequence to every protein, since they should be binned according to it
            let addedSequenceLength =
                let proteinLengthMap =
                    proteinsFromDB
                    |> Array.map (fun (protein,sequence) -> protein, sequence.Length)
                    |> Map.ofArray
                proteins
                |> Array.map (fun ipcis ->
                    ipcis,
                    let groupOfProteins =
                        ipcis.GroupOfProteinIDs
                        |> String.split ';'
                    groupOfProteins
                    |> Array.map (fun protein ->
                        match proteinLengthMap.TryFind protein with
                        | None -> failwith "A protein where you are trying to get the length from isn't present in the database"
                        | Some length -> float length
                    )
                    // lengths are averaged for protein groups
                    |> Array.average
                )
            let combined = Array.append addedSequenceLength proteinsNotMatchedDB
            // sorting treats protein groups as one protein with average length. They are also treated as one protein for total and target counts.
            let sortLength =
                combined
                |> Array.sortBy snd
                |> Array.map fst
            let bins =
                let binSize =
                    ceil (float sortLength.Length / binCount)
                sortLength
                |> Array.chunkBySize (int binSize)
            bins

        // The original paper of Mayu describes a protein as:
        // FP = if only if all its peptides with q <= threshold are decoy
        // TP = at least one of its peptides with q <= threshold is target
        // However, the way mayu estimates it on the program is like this:
        // FP = any protein that contains a decoy psm with q <= threshold
        // TP = any protein that contains a target psm with q <= threshold
        // They do not consider protein containing both decoy and target psms.
        // ProteomIQon currently uses the picked target decoy approach with the following definitions:
        // FP = protein where the score of decoy hits is higher than score of target hits
        // TP = protein where the score of target hits is higher than score of decoy hits
        // Also, mayu estimates q as the empirical (target-decoy) q value.
        // Percolator estimates q as the empirical (target-decoy) q value and adjusted by pi0
        // Mayu extracts the list of TP and FP proteins from PSM level whereas percolator
        // extract the list of TP and FP proteins from peptide level, this avoids redundancy and
        // gives a better calibration since peptide level q values are re-adjusted in percolator.
        // ProteomIQon extracts TP and FP proteins from the result of the picked target decoy approach.
        // This creates sometimes a difference in the number of TP and FP proteins between percolator, Mayu, and ProteomIQon,
        // which causes a slight difference in the estimated protein FDR.

        /// Calculates the expected false positives for every protein bin and sums them up.
        let expectedFP (proteinBin: InferredProteinClassItemScored []) =
            let numberTarget =
                proteinBin
                |> Array.sumBy (fun protein ->
                    match not protein.DecoyHasBetterScore && protein.FoundInDB with
                    | true -> 1.
                    | false -> 0.
                )
            let numberDecoy =
                proteinBin
                 |> Array.sumBy (fun protein ->
                     match protein.DecoyHasBetterScore && protein.FoundInDB with
                     | true -> 1.
                     | false -> 0.
                 )
            let total =
                let notFound =
                    proteinBin
                    |> Array.sumBy (fun protein ->
                        match not protein.FoundInDB with
                        | true -> 1.
                        | false -> 0.
                    )
                notFound + numberTarget + numberDecoy
            // MAYU rounds the number of expected false positives for every bin to the first decimal
            let fpInBin = estimatePi0HG total numberTarget numberDecoy
            fpInBin

    
    /// Calculates the fdr of the data using the MAYU method. The proteinsFromDB is the DB that was used for the inference.
    let calculateFDRwithMAYU (data: InferredProteinClassItemScored []) (proteinsFromDB: (string*string)[]) =
        let proteinBins = MAYU.binProteinsLength data proteinsFromDB 10.
        let estimatedFP =
            proteinBins
            |> Array.fold (fun acc proteinBin -> acc + MAYU.expectedFP proteinBin) 0.
        let targetCount =
            data
            |> Array.sumBy (fun x ->
                match x.DecoyHasBetterScore with
                | true -> 0.
                | false -> 1.
            )
        let fdr =
            if (isNan estimatedFP) || (isInf estimatedFP) || estimatedFP = 0. then
                1.
            elif (estimatedFP / targetCount < 0.) || (estimatedFP / targetCount > 1.) then
                1.
            else
                estimatedFP / targetCount
        fdr

    /// Calculates Decoy/Target ratio
    let calculateFDRwithDecoyTargetRatio (data: InferredProteinClassItemScored []) =
        // Should decoy Hits be doubled?: Target-decoy search strategy for increasedconfidence in large-scale proteinidentifications by mass spectrometry
        let decoyCount  =
            data
            |> Array.sumBy (fun x ->
                match x.DecoyHasBetterScore with
                | true -> 1.
                | false -> 0.
            )
        let targetCount =
            data
            |> Array.sumBy (fun x ->
                match x.DecoyHasBetterScore with
                | true -> 0.
                | false -> 1.
            )
        decoyCount/targetCount


    // Assigns a q value to an InferredProteinClassItemScored
    let assignQValueToIPCIS (qValueF: float -> float) (item: InferredProteinClassItemScored) =
        if item.Decoy then
            createInferredProteinClassItemQValue item (qValueF item.DecoyScore)
        else
            createInferredProteinClassItemQValue item (qValueF item.TargetScore)


<|start_filename|>src/BioFSharp.Mz/XScoring.fs<|end_filename|>
﻿namespace BioFSharp.Mz

open System
open BioFSharp
open BioFSharp.Mz.Peaks
open SearchDB
open Fragmentation
open TheoreticalSpectra
open SearchEngineResult
open FSharpAux.Array
open FSharp.Stats.SpecialFunctions.Factorial

module XScoring =
    
    type MatchingScore = {
        Score   : float 
        N       : int 
        K       : int
        Q       : int
        //Matches : PeakAnnotation list 
        }
    
    let createMatchingScore score n k q = {
     Score = score; N = n ; K = k; Q = q ;//Matches = matches
    }

    [<Struct>]
    type RatedPeak = {
        Peak: float*float
        MoreIntensePeaksWithinWindow: int
        }

    let createRatedPeak peak moreIntensePeaks = {
        Peak = peak; MoreIntensePeaksWithinWindow= moreIntensePeaks }

    let getIntensity rPeak =
        snd rPeak.Peak 
    /// Counts how many peaks present in spectrum data within a a given mzWindow (mzwindow/2 +/- ) have a higher 
    /// intensity then the peak at position spectrumData.[peakIdx]. 
    let private moreIntensePeaksWithin mzWindow (spectrumData:(float*float) []) peakIdx =
        let halfWindowSize = mzWindow / 2. 
        let currentPeakMz,currentPeakIntensity = spectrumData.[peakIdx] 
        let upperBorder = currentPeakMz + halfWindowSize
        let lowerBorder = currentPeakMz - halfWindowSize
        let rightSidePeaks =
            let rec loop currentIdx moreIntensePeaks = 
                if currentIdx = spectrumData.Length then moreIntensePeaks
                else 
                    let tempPeakMz,tempPeakInt = spectrumData.[currentIdx] 
                    if tempPeakMz > upperBorder then moreIntensePeaks
                    elif tempPeakInt > currentPeakIntensity then 
                         loop (currentIdx+1) (moreIntensePeaks+1) 
                    else loop (currentIdx+1) moreIntensePeaks 
            loop peakIdx 0 
        let leftSidePeaks = 
            let rec loop currentIdx moreIntensePeaks = 
                if currentIdx < 0 then moreIntensePeaks
                else 
                    let tempPeakMz,tempPeakInt = spectrumData.[currentIdx] 
                    if tempPeakMz < lowerBorder then moreIntensePeaks
                    elif tempPeakInt > currentPeakIntensity then 
                         loop (currentIdx-1) (moreIntensePeaks+1) 
                    else loop (currentIdx-1) moreIntensePeaks 
            loop peakIdx 0         
        createRatedPeak (currentPeakMz,currentPeakIntensity) (leftSidePeaks+rightSidePeaks) 
    
    /// Applies the function "moreIntensePeaksWithin" to every peak in spectrumData, thus rating every peak by counting how many more abundant peaks
    /// are present within the given mzWindow (mzwindow/2 +/- ). Filters out peaks that are not within the scanLimits and show to have more abundand peaks than 
    /// defined by qMostAbundandpepsMax in their neighbourhood.
    let private ratedSpectrum mzWindow (qMostAbundandPepsMin, qMostAbundandPepsMax) (lowerScanLimit,upperScanLimit) (spectrumData:(float*float) []) = 
        if spectrumData.Length = 0 then [||]
        else
        let filteredAndRanked = 
            spectrumData
            |> Array.mapi   (fun i (mzData,intensityData) -> 
                                moreIntensePeaksWithin mzWindow spectrumData i      
                            )   
            |> Array.filter (fun rPeak -> 
                                let mzData = fst rPeak.Peak  
                                // scanLimits vllt raus
                                mzData >= lowerScanLimit && mzData <= upperScanLimit && rPeak.MoreIntensePeaksWithinWindow < qMostAbundandPepsMax
                            )                        
        filteredAndRanked
          
    /// Converts the fragmentIon ladder to a theoretical spectrum at a given charge state. Filters out all theoretical peaks, that lie
    /// outside the given lower and upperScanLimits.
    let predictOf (lowerScanLimit,upperScanLimit) chargeState (fragments:PeakFamily<TaggedMass.TaggedMass> list)  =
        let predictPeak charge (taggedMass: TaggedMass.TaggedMass) = 
            TaggedPeak.TaggedPeak(taggedMass.Iontype, (Mass.toMZ taggedMass.Mass charge), nan)
        let computePeakFamily charge fragments = 
            let mainPeak = predictPeak charge fragments.MainPeak
            let dependentPeaks =
                fragments.DependentPeaks
                |> List.fold (fun acc dependent -> 
                                 if charge <= 1. then 
                                     predictPeak charge dependent 
                                     :: acc
                                 else 
                                     acc
                              ) []
            createPeakFamily mainPeak dependentPeaks
        let rec recloop ions charge (fragments:PeakFamily<TaggedMass.TaggedMass> list) =
            match fragments with
            | fragments::rest ->  
                if charge <= 1. then 
                    let tempIons = computePeakFamily 1. fragments 
                    recloop (tempIons::ions) charge rest        
                else 

                let tempIons = 
                    [
                    for z = 1 to 2 do 
                        yield computePeakFamily (float z) fragments
                    ]
                recloop (tempIons@ions) charge rest
            | [] -> ions       
        recloop [] chargeState fragments
        //|> List.sortBy (fun peak -> peak.MainPeak.Mz)
        |> List.toArray

    /// Searches for a Peak in ratedSpectrum that matches the target mz (tarmz) within a certain tolerance defined by mzmatchingTolerance.
    let private hasMatchingPeakMZTol (mzMatchingTolerance: float) (tarmz: float) (ratedSpectrum: RatedPeak []) =        
        match binarySearchIndexBy  (fun (rPeak: RatedPeak) ->
                                let mz = fst rPeak.Peak  
                                if ((mz-tarmz) |> abs) <= mzMatchingTolerance then 0 
                                elif tarmz < mz then -1
                                else 1
                            ) ratedSpectrum with
        | x when x >= 0 -> Some ratedSpectrum.[x] 
        | x             ->  None 

    /// Searches for a Peak in ratedSpectrum that matches the target mz (tarmz) within a certain tolerance expressed in parts per million of the tarmz, defined by mzMatchingTolerancePPM.
    let private hasMatchingPeakPPM (mzMatchingTolerancePPM: float) (tarmz: float) (ratedSpectrum: RatedPeak []) = 
        match binarySearchIndexBy (fun (rPeak: RatedPeak) ->
                                let mz = fst rPeak.Peak  
                                if ((mz-tarmz) |> abs) <= Mass.deltaMassByPpm mzMatchingTolerancePPM tarmz then 0 
                                elif tarmz < mz then -1
                                else 1
                            ) ratedSpectrum with
        | x when x >= 0 -> Some ratedSpectrum.[x] 
        | _             -> None 


    [<Struct>]
    type CountedMatches = {
        Q                       : int
        N_BIons                 : int
        N_YIons                 : int
        N_Neutral               : int
        K_BIons                 : int
        K_YIons                 : int
        K_Neutral               : int
        MatchedSum              : float
        }

    let createCountedMatches q n_BIons n_YIons n_Neutral k_BIons k_YIons k_Neutral matchedSum = {
        Q           = q                       
        N_BIons     = n_BIons                 
        N_YIons     = n_YIons                 
        N_Neutral   = n_Neutral               
        K_BIons     = k_BIons                 
        K_YIons     = k_YIons                 
        K_Neutral   = k_Neutral               
        MatchedSum  = matchedSum  
        }

    let findMaxIterator qMostAbundandPepsMin countArrayL (rPeak:RatedPeak) = 
        match (1+rPeak.MoreIntensePeaksWithinWindow)-qMostAbundandPepsMin with 
        | x when x <= 0 -> 0
        | x when x >= countArrayL-1 -> x
        | x -> x      

    let raiseCountsAndHits qMostAbundandPepsMin (countArray: CountedMatches []) (theoPeak:TaggedPeak.TaggedPeak) (rPeak:RatedPeak) =
        let maxIterator = findMaxIterator qMostAbundandPepsMin countArray.Length rPeak
        match theoPeak.Iontype with
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.B) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                if i >= maxIterator then
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                        (currentMatchingCount.N_BIons+1)                 
                                        (currentMatchingCount.N_YIons)                 
                                        (currentMatchingCount.N_Neutral)               
                                        (currentMatchingCount.K_BIons+1)                 
                                        (currentMatchingCount.K_YIons)                 
                                        (currentMatchingCount.K_Neutral)               
                                        (currentMatchingCount.MatchedSum+getIntensity rPeak)
                else 
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                       (currentMatchingCount.N_BIons+1)                 
                                       (currentMatchingCount.N_YIons)                 
                                       (currentMatchingCount.N_Neutral)               
                                       (currentMatchingCount.K_BIons)                 
                                       (currentMatchingCount.K_YIons)                 
                                       (currentMatchingCount.K_Neutral)               
                                       (currentMatchingCount.MatchedSum) 
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.Y) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                if i >= maxIterator then
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                        (currentMatchingCount.N_BIons)                 
                                        (currentMatchingCount.N_YIons+1)                 
                                        (currentMatchingCount.N_Neutral)               
                                        (currentMatchingCount.K_BIons)                 
                                        (currentMatchingCount.K_YIons+1)                 
                                        (currentMatchingCount.K_Neutral)               
                                        (currentMatchingCount.MatchedSum+getIntensity rPeak)
                else 
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                       (currentMatchingCount.N_BIons)                 
                                       (currentMatchingCount.N_YIons+1)                 
                                       (currentMatchingCount.N_Neutral)               
                                       (currentMatchingCount.K_BIons)                 
                                       (currentMatchingCount.K_YIons)                 
                                       (currentMatchingCount.K_Neutral)               
                                       (currentMatchingCount.MatchedSum) 
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.Neutral) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                if i >= maxIterator then
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                        (currentMatchingCount.N_BIons)                 
                                        (currentMatchingCount.N_YIons)                 
                                        (currentMatchingCount.N_Neutral+1)               
                                        (currentMatchingCount.K_BIons)                 
                                        (currentMatchingCount.K_YIons)                 
                                        (currentMatchingCount.K_Neutral+1)               
                                        (currentMatchingCount.MatchedSum+getIntensity rPeak)   
                else 
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                       (currentMatchingCount.N_BIons)                 
                                       (currentMatchingCount.N_YIons)                 
                                       (currentMatchingCount.N_Neutral+1)               
                                       (currentMatchingCount.K_BIons)                 
                                       (currentMatchingCount.K_YIons)                 
                                       (currentMatchingCount.K_Neutral)               
                                       (currentMatchingCount.MatchedSum)    
        | _ -> 
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                if i >= maxIterator then
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                        (currentMatchingCount.N_BIons)                 
                                        (currentMatchingCount.N_YIons)                 
                                        (currentMatchingCount.N_Neutral)               
                                        (currentMatchingCount.K_BIons)                 
                                        (currentMatchingCount.K_YIons)                 
                                        (currentMatchingCount.K_Neutral)               
                                        (currentMatchingCount.MatchedSum)                
                else 
                    countArray.[i] <- createCountedMatches 
                                        currentMatchingCount.Q                       
                                       (currentMatchingCount.N_BIons)                 
                                       (currentMatchingCount.N_YIons)                 
                                       (currentMatchingCount.N_Neutral)               
                                       (currentMatchingCount.K_BIons)                 
                                       (currentMatchingCount.K_YIons)                 
                                       (currentMatchingCount.K_Neutral)               
                                       (currentMatchingCount.MatchedSum) 

    let raiseOnlyCounts (countArray: CountedMatches []) (theoPeak:TaggedPeak.TaggedPeak) =
        match theoPeak.Iontype with
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.B) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                countArray.[i] <- createCountedMatches 
                                    currentMatchingCount.Q                       
                                    (currentMatchingCount.N_BIons+1)                 
                                    (currentMatchingCount.N_YIons)                 
                                    (currentMatchingCount.N_Neutral)               
                                    (currentMatchingCount.K_BIons)                 
                                    (currentMatchingCount.K_YIons)                 
                                    (currentMatchingCount.K_Neutral)               
                                    (currentMatchingCount.MatchedSum) 
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.Y) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                countArray.[i] <- createCountedMatches 
                                    currentMatchingCount.Q                       
                                    (currentMatchingCount.N_BIons)                 
                                    (currentMatchingCount.N_YIons+1)                 
                                    (currentMatchingCount.N_Neutral)               
                                    (currentMatchingCount.K_BIons)                 
                                    (currentMatchingCount.K_YIons)                 
                                    (currentMatchingCount.K_Neutral)               
                                    (currentMatchingCount.MatchedSum) 
        | iontype when iontype.HasFlag(Ions.IonTypeFlag.Neutral) ->
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                countArray.[i] <- createCountedMatches 
                                    currentMatchingCount.Q                       
                                    (currentMatchingCount.N_BIons)                 
                                    (currentMatchingCount.N_YIons)                 
                                    (currentMatchingCount.N_Neutral+1)               
                                    (currentMatchingCount.K_BIons)                 
                                    (currentMatchingCount.K_YIons)                 
                                    (currentMatchingCount.K_Neutral)               
                                    (currentMatchingCount.MatchedSum)   
        | _ -> 
            for i = 0 to countArray.Length-1 do
                let currentMatchingCount = countArray.[i] 
                countArray.[i] <- createCountedMatches 
                                    currentMatchingCount.Q                       
                                    (currentMatchingCount.N_BIons)                 
                                    (currentMatchingCount.N_YIons)                 
                                    (currentMatchingCount.N_Neutral)               
                                    (currentMatchingCount.K_BIons)                 
                                    (currentMatchingCount.K_YIons)                 
                                    (currentMatchingCount.K_Neutral)               
                                    (currentMatchingCount.MatchedSum)  
        
    /// Iterates theoSpect and counts the number of common peaks in ratedspectrum. A peak is considered as common if the mz value lies within a user defined matching tolerance. This procedure
    /// is carried out (qMostAbundandPepsMax-qMostAbundandPepsMin+1) times every time with a different qThreshold to account for peak intensity. A peak is counted as matching if it has
    /// less more abundant peaks in the neighbourhood than the current qThreshold.
    let private countMatches (lowerScanLimit,upperScanLimit) (qMostAbundandPepsMin, qMostAbundandPepsMax) (mzMatchingTolerance: float) (theoSpect: PeakFamily<TaggedPeak.TaggedPeak> []) (ratedSpectrum: RatedPeak []) =    
        let countArray = Array.init (qMostAbundandPepsMax-qMostAbundandPepsMin+1) (fun i -> createCountedMatches (qMostAbundandPepsMin+i) 0 0 0 0 0 0 0.)   
        let rec findMatches (mzMatchingTolerance: float) (theoSpect: PeakFamily<TaggedPeak.TaggedPeak> []) (ratedSpectrum: RatedPeak []) idx =
            if idx = theoSpect.Length then
                countArray
            else
                let currentPeakFam = theoSpect.[idx]
                if currentPeakFam.MainPeak.Mz >= lowerScanLimit && currentPeakFam.MainPeak.Mz <= upperScanLimit then                
                    match hasMatchingPeakPPM mzMatchingTolerance currentPeakFam.MainPeak.Mz ratedSpectrum with 
                    | Some rPeak ->
                        raiseCountsAndHits qMostAbundandPepsMin countArray currentPeakFam.MainPeak rPeak
                        currentPeakFam.DependentPeaks 
                        |> List.iter ( fun dependentAnnotPeak ->
                                        // Peak not within bounds loop further. 
                                        if dependentAnnotPeak.Mz >= lowerScanLimit && dependentAnnotPeak.Mz <= upperScanLimit then 
                                             match hasMatchingPeakPPM mzMatchingTolerance dependentAnnotPeak.Mz ratedSpectrum with 
                                             //Successful match: raise n and k
                                             | Some rPeak -> raiseCountsAndHits qMostAbundandPepsMin countArray dependentAnnotPeak rPeak
                                             | None       -> raiseOnlyCounts countArray dependentAnnotPeak 
                                        )  
                        findMatches mzMatchingTolerance theoSpect ratedSpectrum (idx+1) 
                    | None ->
                        raiseOnlyCounts countArray currentPeakFam.MainPeak
                        findMatches mzMatchingTolerance theoSpect ratedSpectrum (idx+1)
                else 
                    // Peak not within bounds 
                    findMatches mzMatchingTolerance theoSpect ratedSpectrum (idx+1)
        findMatches mzMatchingTolerance theoSpect ratedSpectrum 0                          

    /// Helper function to calculate the AndromedaScore
    let private lnProb n k lnp lnq = 
        let s1 = -k * lnp - (n - k) * lnq - FSharp.Stats.SpecialFunctions.Gamma.gammaLn(n + 1.) 
                    + FSharp.Stats.SpecialFunctions.Gamma.gammaLn(k + 1.) + FSharp.Stats.SpecialFunctions.Gamma.gammaLn(n - k + 1.) 
        s1 
    
    /// Calculates the andromedaLike Score based on offered peaks n and matched peaks k, as well as the q value threshold topx
    let private scoreFuncImpl n k topx = 
        let fTopx = float topx 
        let p1 = Math.Min(fTopx/100.0,0.5)
        let lnp = Math.Log(p1)
        let lnq = Math.Log(1.0-p1)
        let mutable acc = 0.
        for i = k to n do 
            acc <- acc + Math.Exp(-lnProb (float n) (float i) lnp lnq)
        acc <- -Math.Log(acc)
        10. * acc / Math.Log(10.)

    /// Correction factor applided to the andromeda score based on the observed precursor Mz.
    /// Values are taken from the original andromeda release.
    let private mzCorrection mass = 
        0.024*(mass - 600.)

    /// Correction factor applided to the andromeda score based on the amount modifications.
    /// Values are taken from the original andromeda release.
    let private modCorrection nMod = 
        match nMod with 
        | 0  -> 42.
        | 1  -> 28.
        | 2  -> 22.
        | 3  -> 16.
        | 4  -> 9.
        | 5  -> 5.
        | 6  -> 2.
        | _  -> 0.
    
    /// Correction factor applided to the andromeda score based on the amount of missCleavages.
    /// Values are taken from the original andromeda release.
    let private cleavageCorrection nCleave consecutiveCleavage =
        match nCleave with 
        | 0 -> 53.2
        | 1 -> 
            if consecutiveCleavage then 
                42.1
            else 
                31.1
        | 2 -> 17. 
        |_  -> 0.

    let optimizeAndromedaScore precursorMZ countedMatches = 
        let massCorr = mzCorrection precursorMZ //(Mass.toMZ lookUpResult.Mass charge)
        let modCorr  = modCorrection 0
        let cleavageCorr = cleavageCorrection 0 false
        let nWithoutNeutral = countedMatches.N_BIons+countedMatches.N_YIons
        let kWithoutNeutral = countedMatches.K_BIons+countedMatches.K_YIons
        let nWithNeutral = countedMatches.N_BIons+countedMatches.N_YIons+countedMatches.N_Neutral
        let kWithNeutral = countedMatches.K_BIons+countedMatches.K_YIons+countedMatches.K_Neutral
        let rawScore1 = scoreFuncImpl nWithoutNeutral kWithoutNeutral countedMatches.Q
        let rawScore2 = scoreFuncImpl nWithNeutral kWithNeutral countedMatches.Q
        if rawScore1 > rawScore2 then 
            let finalScore = 
                let tmp = (rawScore1 + massCorr + modCorr + cleavageCorr - 100.)
                if tmp < 0. then 0. else tmp 
            createMatchingScore finalScore nWithoutNeutral kWithoutNeutral countedMatches.Q
        else 
            let finalScore = 
                let tmp = (rawScore2 + massCorr + modCorr + cleavageCorr - 100.)
                if tmp < 0. then 0. else tmp 
            createMatchingScore finalScore nWithNeutral kWithNeutral countedMatches.Q 
    
    /// Calculates the X!TandemHyperScore
    let calcHyperScore (countedMatches:CountedMatches) =
        let res = log countedMatches.MatchedSum + factorialLn countedMatches.K_BIons + factorialLn countedMatches.K_YIons + factorialLn countedMatches.K_Neutral
        Math.Max(0.,res)

    /// Computes the andromedalike score of a theoretical spectrum vs. the ratedSpectrum. Subsequently applies mzCorrection - modCorrection and cleavageCorrection to the score. 
    let private scoreTheoVsRecordedSpec (lowerScanLimit,upperScanLimit) (qMostAbundandPepsMin, qMostAbundandPepsMax)  matchingTolPPM  precursorMZ (*charge*) theoSpec (ratedSpectrum: RatedPeak []) =     
        let matches = countMatches (lowerScanLimit,upperScanLimit) (qMostAbundandPepsMin, qMostAbundandPepsMax) matchingTolPPM theoSpec ratedSpectrum
        let androScore = 
            matches
            |> Array.map (optimizeAndromedaScore precursorMZ)
            |> Array.maxBy (fun x -> x.Score) 
        let bestXTandem = 
            matches
            |> Array.maxBy (fun x -> x.Q) 
        let xTandemScore = 
            bestXTandem
            |> calcHyperScore
        let matchSum = bestXTandem.MatchedSum
        androScore,xTandemScore,matchSum

    /// Converts the fragment ion ladders to a theoretical Sequestlike spectrum at a given charge state. 
    /// Subsequently, the spectrum is binned to the nearest mz bin (binwidth = 1 Da). Filters out peaks 
    /// that are not within the scanLimits.
    let getTheoSpecs scanlimits chargeState (possiblePeptideInfos:list<LookUpResult<AminoAcids.AminoAcid>*FragmentMasses>) =
        TheoreticalSpectra.getTheoSpecs predictOf scanlimits chargeState possiblePeptideInfos
 
    /// Calculates the AndromedaLike scores for all theoretical spectra.
    let calcAndromedaAndXTandemScore qMinAndMax scanlimits matchingTolPPM (spectrum:PeakArray<_>) scanTime chargeState isolationWindowTargetMz (theoreticalSpectra:list<TheoreticalSpectrum<PeakFamily<TaggedPeak.TaggedPeak>[]>> ) spectrumID =
        
        let fCharge = float chargeState
        //
        let binnedRecSpec = 
            spectrum 
            |> Array.map (fun peak -> peak.Mz,peak.Intensity)
            |> ratedSpectrum 100. qMinAndMax scanlimits
        //
        let ms_mass = Mass.ofMZ isolationWindowTargetMz fCharge
        let androScores,xTandemScores = 
            theoreticalSpectra 
            |> List.fold (fun (andro,xTandem) theoreticalSpectrum -> 
                                let lookUpResult = theoreticalSpectrum.LookUpResult
                                let seqL = lookUpResult.BioSequence.Length
 
                                // Calculates score for target sequence
                                let theoSpec = theoreticalSpectrum.TheoSpec
                                let targetScoreAndro,targetScoreXTandem,matchSum = scoreTheoVsRecordedSpec scanlimits qMinAndMax matchingTolPPM isolationWindowTargetMz theoSpec binnedRecSpec
                                let targetResultAndro = createSearchEngineResult SearchEngineResult.SearchEngine.AndromedaLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod true scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass seqL targetScoreAndro.Score nan nan
                                let targetResultXTandem = createSearchEngineResult SearchEngineResult.SearchEngine.XTandemLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod true scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass seqL targetScoreXTandem nan nan
                                
                                let theoSpec_Decoy = theoreticalSpectrum.DecoyTheoSpec
                                let decoyScoreAndro,decoyScoreXTandem,matchSum = scoreTheoVsRecordedSpec scanlimits qMinAndMax matchingTolPPM isolationWindowTargetMz theoSpec_Decoy binnedRecSpec
                                let decoyResultAndro = createSearchEngineResult SearchEngineResult.SearchEngine.AndromedaLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod false scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass seqL decoyScoreAndro.Score nan nan
                                let decoyResultXTandem = createSearchEngineResult SearchEngineResult.SearchEngine.XTandemLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod false scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass seqL decoyScoreXTandem nan nan
                                
                                targetResultAndro :: decoyResultAndro :: andro,targetResultXTandem :: decoyResultXTandem :: xTandem      
                     ) ([],[])
        let androScoresFinal =
            calcNormDeltaBestToRest (androScores |> List.sortBy (fun sls -> - sls.Score))       
            |> calcNormDeltaNext
        let xTandemScoresFinal =
            calcNormDeltaBestToRest (xTandemScores |> List.sortBy (fun sls -> - sls.Score))       
            |> calcNormDeltaNext
        androScoresFinal,xTandemScoresFinal 
   
   

<|start_filename|>src/BioFSharp.Mz/FDRControl.fs<|end_filename|>
﻿namespace BioFSharp.Mz

open FSharp.Stats
open FSharpAux
open FSharp.Stats.Fitting
open FSharp.Stats.Fitting.NonLinearRegression

module FDRControl = 

    module MAYU =

        // FDR estimation using MAYU
        // Code form 'stirlingLogFactorial' to 'estimatePi0HG' translated from percolator 'ProteinFDREstimator.cpp'
        let private stirlingLogFacorial (n: float) =
            log(sqrt(2. * pi  *n)) + n * log(n) - n

        let private exactLogFactorial (n: float) =
            let rec loop i log_fact =
                if i > n then
                    log_fact
                else
                    let new_log_fact = log_fact + (log i)
                    loop (i + 1.) new_log_fact
            loop 2. 0.

        let private logFactorial (n: float) =
            if n < 1000. then
                exactLogFactorial n
            else
                stirlingLogFacorial n

        let private logBinomial (n: float) (k: float) =
            (logFactorial n) - (logFactorial k) - (logFactorial (n - k))

        let private hypergeometric (x: float) (n: float) (w: float) (d: float) =
            //natural logarithm of the probability
            if (d > 0.) then
                exp((logBinomial w x) + (logBinomial (n - w) (d - x)) - (logBinomial n d))
            else 0.

        /// Estimates the false positives given the total number of entries, the number of target hits and the number of decoy hits
        let estimatePi0HG (n: float) (targets: float) (cf: float) =
            let logprob = 
                [0. .. cf]
                |> List.fold (fun acc fp ->
                    let tp = targets - fp
                    let w = n - tp
                    let prob = hypergeometric fp n w cf
                    prob::acc
                ) []
                |> List.rev
            let sum = logprob |> List.sum
            let logprob_Norm =
                logprob
                |> List.map (fun x -> x / sum)
            // MAYU rounds here to first decimal
            let expectation_value_FP_PID =
                logprob_Norm
                |> List.foldi (fun i acc x -> acc + x * (float i)) 0.
            if (isNan expectation_value_FP_PID) || (isInf expectation_value_FP_PID) then
                0.
            else
                expectation_value_FP_PID

    type ScoreTargetDecoyCount =
        {
            Score      : float
            DecoyCount : float
            TargetCount: float
        }
     
    /// Gives the decoy and target count at a specific score
    let createScoreTargetDecoyCount score decoyCount targetCount =
        {
            Score       = score
            DecoyCount  = decoyCount
            TargetCount = targetCount
        }
    
    /// returns scores, pep, q
    let binningFunction bandwidth pi0 (scoreF: 'A -> float) (isDecoyF: 'A -> bool) (data:'A[])  = 
        let totalDecoyProportion = 
            let decoyCount = 
                Array.filter isDecoyF data
                |> Array.length
                |> float
            let totalCount =
                data
                |> Array.length
                |> float
            1. / (2. * decoyCount / totalCount)
        data
        |> Array.groupBy (fun s -> floor (scoreF s / bandwidth))
        |> Array.sortBy fst
        |> Array.map (fun (k,values)->
            let median     = values |> Array.map scoreF |> Array.average
            let totalCount = values |> Array.length |> float
            let decoyCount = values |> Array.filter isDecoyF |> Array.length |> float |> (*) totalDecoyProportion
            median,totalCount,decoyCount
        )
        |> fun a ->
            a
            |> Array.mapi (fun i (median,totalCountBin,decoyCountBin) ->
                            /// TODO: Accumulate totalCount + totalDecoyCount beforeHand and skip the time intensive mapping accross the array in each iteration.
                            let _,totalCountRight,decoyCountRight = a.[i..a.Length-1] |> Array.reduce (fun (x,y,z) (x',y',z') -> x+x',y+y',z+z')
                            (median,(pi0 * 2. * decoyCountBin / totalCountBin),(pi0 * 2. * decoyCountRight / totalCountRight))
                          )
        |> Array.sortBy (fun (score,pep,q) -> score) 
        |> Array.unzip3 
        |> fun (score,pep,q) -> vector score, vector pep, vector q
    


    /// Input for QValue calulation
    type QValueInput =
        {
            Score    : float
            IsDecoy  : bool
        }

    let createQValueInput score isDecoy =
        {
            Score     = score
            IsDecoy   = isDecoy
        }

    /// Gives a function to calculate the q value for a score in a dataset using Lukas method and Levenberg Marguardt fitting
    let calculateQValueLogReg fdrEstimate bandwidth (data: 'a []) (isDecoy: 'a -> bool) (decoyScoreF: 'a -> float) (targetScoreF: 'a -> float) =
        // Input for q value calculation
        let createTargetDecoyInput =
            data
            |> Array.map (fun item ->
                if isDecoy item then
                    createQValueInput (decoyScoreF item) true
                else
                    createQValueInput (targetScoreF item) false
            )

        let scores,pep,qVal =
            binningFunction bandwidth fdrEstimate (fun (x: QValueInput) -> x.Score) (fun (x: QValueInput) -> x.IsDecoy) createTargetDecoyInput
            |> fun (scores,pep,qVal) -> scores.ToArray(), pep.ToArray(), qVal.ToArray()

        //Chart.Point (scores,qVal)
        //|> Chart.Show

        // gives a range of 1 to 30 for the steepness. This can be adjusted depending on the data, but normally it should lie in this range
        let initialGuess =
            Fitting.NonLinearRegression.LevenbergMarquardtConstrained.initialParamsOverRange scores qVal [|1. .. 30.|]
            |> Array.map (fun guess -> Table.lineSolverOptions guess)

        // performs Levenberg Marguardt Constrained algorithm on the data for every given initial estimate with different steepnesses and selects the one with the lowest RSS
        let estimate =
            initialGuess
            |> Array.map (fun initial ->
                if initial.InitialParamGuess.Length > 3 then failwith "Invalid initial param guess for Logistic Function"
                let lowerBound =
                    initial.InitialParamGuess
                    |> Array.map (fun param -> param - (abs param) * 0.1)
                    |> vector
                let upperBound =
                    initial.InitialParamGuess
                    |> Array.map (fun param -> param + (abs param) * 0.1)
                    |> vector
                LevenbergMarquardtConstrained.estimatedParamsWithRSS Table.LogisticFunctionDescending initial 0.001 10.0 lowerBound upperBound scores qVal
            )
            |> Array.filter (fun (param,rss) -> not (param |> Vector.exists System.Double.IsNaN))
            |> Array.minBy snd
            |> fst

        let logisticFunction = Table.LogisticFunctionDescending.GetFunctionValue estimate
        logisticFunction



    /// Gives a function to calculate the q value for a score in a dataset using Storeys method
    let calculateQValueStorey (data: 'a[]) (isDecoy: 'a -> bool) (decoyScoreF: 'a -> float) (targetScoreF: 'a -> float) =
        // Gives an array of scores with the frequency of decoy and target hits at that score
        let scoreFrequencies =
            data
            |> Array.map (fun x ->
                if isDecoy x then
                    decoyScoreF x, true
                else
                    targetScoreF x, false
            )
            // groups by score
            |> Array.groupBy fst
            // counts occurences of targets and decoys at that score
            |> Array.map (fun (score,scoreDecoyInfo) ->
                let decoyCount =
                    scoreDecoyInfo
                    |> Array.sumBy (fun (score, decoyInfo) ->
                        match decoyInfo with
                        | true -> 1.
                        | false -> 0.
                    )
                let targetCount =
                    scoreDecoyInfo
                    |> Array.sumBy (fun (score, decoyInfo) ->
                        match decoyInfo with
                        | true -> 0.
                        | false -> 1.
                    )
                createScoreTargetDecoyCount score decoyCount targetCount
            )
            |> Array.sortByDescending (fun x -> x.Score)

        // Goes through the list and assigns each protein a "q value" by dividing total decoy hits so far through total target hits so far
        let reverseQVal =
            scoreFrequencies
            |> Array.fold (fun (acc: (float*float*float*float) list) scoreCounts ->
                let _,_,decoyCount,targetCount = acc.Head
                // Should decoy hits be doubled?
                // accumulates decoy hits
                let newDecoyCount  = decoyCount + scoreCounts.DecoyCount(* * 2.*)
                // accumulates target hits
                let newTargetCount = targetCount + scoreCounts.TargetCount
                let newQVal =
                    let nominator =
                        if newTargetCount > 0. then newTargetCount
                        else 1.
                    newDecoyCount / nominator
                (scoreCounts.Score, newQVal, newDecoyCount, newTargetCount):: acc
            ) [0., 0., 0., 0.]
            // removes last part of the list which was the "empty" initial entry
            |> fun list -> list.[.. list.Length-2]
            |> List.map (fun (score, qVal, decoyC, targetC) -> score, qVal)

        //Assures monotonicity by going through the list from the bottom to top and assigning the previous q value if it is smaller than the current one
        let score, monotoneQVal =
            if reverseQVal.IsEmpty then
                failwith "Reverse qvalues in Storey calculation are empty"
            let head::tail = reverseQVal
            tail
            |> List.fold (fun (acc: (float*float) list) (score, newQValue) ->
                let _,qValue = acc.Head
                if newQValue > qValue then
                    (score, qValue)::acc
                else
                    (score, newQValue)::acc
            )[head]
            |> Array.ofList
            |> Array.sortBy fst
            |> Array.unzip
        // Linear Interpolation
        let linearSplineCoeff = Interpolation.LinearSpline.initInterpolateSorted score monotoneQVal
        // takes a score from the dataset and assigns it a q value
        let interpolation = Interpolation.LinearSpline.interpolate linearSplineCoeff
        interpolation

<|start_filename|>src/BioFSharp.Mz/SequestLike.fs<|end_filename|>
﻿namespace BioFSharp.Mz

open BioFSharp
open SearchDB
open Fragmentation
open TheoreticalSpectra
open FSharpAux
open SearchEngineResult
open Ions
open FSharp.Stats
//open MathNet.Numerics
//open MathNet.Numerics.LinearAlgebra.Double

module SequestLike =

    /// normalize the intensities within a window to maximum of the window
    /// Attention shortens the array  (cuts)
    let windowNormalizeIntensities (intensities:Vector<float>) (numberOfWindows:int) =
        // finds max within range
        let rec findMax (array:Vector<float>) (cMax:float) (lowerLimit:int) (counter:int)  =
            if counter < lowerLimit then
                cMax
            else
                let nMax = max array.[counter] cMax
                findMax array nMax lowerLimit (counter - 1)
        // sqrt intensities and normalize to maximum within lower-upper limit
        let normToMaxSqrtInPlace (array:Vector<float>) (lowerLimit:int) (upperLimit:int)  =
            let nMax = sqrt (findMax array 0.0 lowerLimit upperLimit)
            for i = lowerLimit to upperLimit do
                if nMax > 0. then                
                    array.[i] <- (sqrt array.[i]) / nMax 
                else
                    array.[i] <- 0.
        let windowSize =  (intensities.Length / numberOfWindows)    
        let tmpIntensities =            
            Vector.init (windowSize*numberOfWindows) (fun i -> intensities.[i])            
        for i = 1 to numberOfWindows do
            //printfn "window: %i lower: %i counter: %i " i (windowSize * (i - 1)) (windowSize * i - 1)
            normToMaxSqrtInPlace tmpIntensities (windowSize * (i - 1)) (windowSize * i - 1) 
        tmpIntensities
    
    /// Predicts the intensity of a theoretical peak based on the given charge and iontype.
    let predictIntensitySimpleModel (ionType:Ions.IonTypeFlag) (charge:float) =
        match ionType with 
        | ionT when hasFlag Ions.IonTypeFlag.Unknown    ionT -> 1.  / charge
        | ionT when hasFlag Ions.IonTypeFlag.Diagnostic ionT -> 1.  / charge
        | ionT when hasFlag Ions.IonTypeFlag.Neutral    ionT -> 0.6 / charge
        | ionT when hasFlag Ions.IonTypeFlag.Immonium   ionT -> 0.6 / charge
        | ionT when hasFlag Ions.IonTypeFlag.lossNH3    ionT -> 0.2 / charge
        | ionT when hasFlag Ions.IonTypeFlag.lossH2O    ionT -> 0.2 / charge
        | ionT when hasFlag Ions.IonTypeFlag.A          ionT -> 0.2 / charge   
        | ionT when hasFlag Ions.IonTypeFlag.B          ionT -> 1.  / charge
        | ionT when hasFlag Ions.IonTypeFlag.C          ionT -> 0.2 / charge
        | ionT when hasFlag Ions.IonTypeFlag.X          ionT -> 0.2 / charge
        | ionT when hasFlag Ions.IonTypeFlag.Y          ionT -> 1.  / charge
        | ionT when hasFlag Ions.IonTypeFlag.Z          ionT -> 0.2 / charge
        | _                                                  -> 0.2 / charge

    let setVectorInplace (vector:Vector<float>) charge maxIndex lowerScanLimit (taggedMass:TaggedMass.TaggedMass) = 
        let index = int(System.Math.Round (Mass.toMZ taggedMass.Mass charge) ) - lowerScanLimit 
        if index < maxIndex-1 && index > -1 then
            vector.[index] <- max vector.[index] (predictIntensitySimpleModel (taggedMass.Iontype) charge) 

    /// Converts the fragment ion ladder to a theoretical Sequestlike spectrum at a given charge state. 
    let theoSpecOf (lowerScanLimit,upperScanLimit) (maxcharge:float) (fragments:PeakFamily<TaggedMass.TaggedMass> list) =
        let lowerScanLimit = int lowerScanLimit
        let upperScanLimit = int upperScanLimit
        let maxIndex = upperScanLimit - lowerScanLimit + 1        
        let vector = Vector.create (maxIndex-1) 0.
        fragments 
        |> List.iter (fun p ->  
            let peaks = p.MainPeak::p.DependentPeaks
            [1. .. maxcharge]  
            |> List.iter (fun ch -> peaks |> List.iter (setVectorInplace vector ch maxIndex lowerScanLimit) ) 
            )
        vector

    /// Computes the autocorrelation of the vector +/- plusMinusMaxDelay.
    let autoCorrelation (plusMinusMaxDelay:int) (vector:Vector<float>) =
        let shifted (vector:Vector<float>) (tau:int) =
            vector
            |> Vector.mapi
                (fun i x ->
                    let index = i - tau
                    if (index < 0) || (index > vector.Length - 1) then 
                        0.
                    else
                        vector.[index] )
        let rec accumVector (accum) (state:int) (max:int) =
            if state = max then
                accum
            else
                accumVector (accum + (shifted vector state)) (state - 1) (max)
        let emtyVector = Vector.zero vector.Length
        let plus  = accumVector emtyVector (plusMinusMaxDelay) (1)
        let minus = accumVector emtyVector (-1) (-plusMinusMaxDelay)
        (plus + minus)
        |> Vector.map (fun x ->  x / (float plusMinusMaxDelay * 2.))
       
    /////
    //let private createShiftedMatrix (plusMinusMaxDelay:int) (array:float[]) =
    //    let colNumber = array.Length
    //    Array2D.init (plusMinusMaxDelay * 2 + 1) colNumber
    //        (fun i ii ->
    //            let ni = (i - plusMinusMaxDelay)
    //            let index = ii - ni
    //            if (index < 0) || (index > colNumber - 1) then 
    //                0.
    //            else
    //                array.[index] )

    /// Converts the fragment ion ladder to a theoretical Sequestlike spectrum at a given charge state. 
    /// Subsequently, the spectrum is binned to the nearest mz bin (binwidth = 1 Da). Filters out peaks 
    /// that are not within the scanLimits.
    let predictOf (lowerScanLimit,upperScanLimit) charge ionSeries =         
        theoSpecOf (lowerScanLimit,upperScanLimit) charge ionSeries      
        

//    /// Measured spectrum to sequest-like normalized intensity array
//    /// ! Uses 10 as number of windows for window normalization (like in original sequest algorithm)    
//    let spectrumToNormalizedIntensityArray (scanlimits:int*int) (spectrum:PeakArray<_>) =
//        let lowerScanLimit,upperScanLimit = scanlimits    
//        let si  = PeakArray.peaksToNearestUnitDaltonBin spectrum lowerScanLimit upperScanLimit
//        let nsi = windowNormalizeIntensities si 10 // |> Seq.toArray
//        nsi

    /// Measured spectrum to sequest-like normalized intensity array
    /// minus auto-correlation (delay 75 -> like in original sequest algorithm)
    /// ! Uses 10 as number of windows for window normalization (like in original sequest algorithm)    
    let spectrumToIntensityArrayMinusAutoCorrelation (lowerScanLimit,upperScanLimit) (spectrum:PeakArray<_>) =
        let si  = PeakArray.peaksToNearestUnitDaltonBinVector spectrum lowerScanLimit upperScanLimit
        let nsi = windowNormalizeIntensities si 10    
        let nsi' = autoCorrelation 75 nsi
        (nsi - nsi')  


    /// Calculates the Cross-Correlation of p_nis and ms_nis
    let calcXCorr (p_nis:Vector<float>) (ms_nis:Vector<float>) =
        let tmp = Vector.dot p_nis ms_nis 
        if tmp < 0. then 0. else tmp 

    /// Substracts Auto Correlation of measured normalized intensity spectrum (ms_nis). Subsequently, 
    /// the Cross-Correlation of the predicted normalized intensity spectrum (p_nis) and ms_nis is calculated. 
    let scoreSingle (p_nis:Vector<float>) (ms_nis:Vector<float>)  =
        let ms_nis' = autoCorrelation 75 ms_nis
        (ms_nis - ms_nis') 
        |> calcXCorr p_nis 
        
    /// Converts the fragment ion ladders to a theoretical Sequestlike spectrum at a given charge state. 
    /// Subsequently, the spectrum is binned to the nearest mz bin (binwidth = 1 Da). Filters out peaks 
    /// that are not within the scanLimits.
    let getTheoSpecs scanlimits chargeState (possiblePeptideInfos:list<LookUpResult<AminoAcids.AminoAcid>*FragmentMasses>) =
        TheoreticalSpectra.getTheoSpecs predictOf scanlimits chargeState possiblePeptideInfos

    /// Calculates the SequestLike Scores for all peptides (possiblePeptideInfos).
    let calcSequestLikeScoresRevDecoy calcIonSeries (massfunction:Formula.Formula -> float) (scanlimits) (spectrum:PeakArray<_>) scanTime chargeState isolationWindowTargetMz (possiblePeptideInfos:list<LookUpResult<AminoAcids.AminoAcid>>) spectrumID =                             
        // measured normailzed intensity array (spectrum) minus auto-correlation
        let ms_nis =  spectrumToIntensityArrayMinusAutoCorrelation scanlimits spectrum
        // float charge
        let fCharge = float chargeState
        // measured mass
        let ms_mass = Mass.ofMZ isolationWindowTargetMz fCharge
        //
        let ides = 
            possiblePeptideInfos 
            |> List.fold (fun acc lookUpResult -> 
                              let sequence = lookUpResult.BioSequence        
                              let ionSeries = calcIonSeries massfunction sequence    
                              //predicted  normalized intensity array (spectrum) 
                              let p_nis = predictOf scanlimits fCharge ionSeries 
                              let xcorr = calcXCorr p_nis ms_nis
                              let targetScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod true scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass  (List.length sequence) xcorr nan nan
                              
                              let revPeptide_decoy = sequence |> List.rev
                              let ionSeries_decoy  = calcIonSeries massfunction revPeptide_decoy
                              let p_nis_decoy      = predictOf scanlimits fCharge ionSeries_decoy 
                              let xcorr_decoy      = calcXCorr p_nis_decoy ms_nis
                              let decoyScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod false scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass revPeptide_decoy.Length xcorr_decoy nan nan
                              targetScore::decoyScore::acc  
                     ) []
        calcNormDeltaBestToRest (ides |> List.sortBy (fun sls -> - sls.Score))        
        |> calcNormDeltaNext

    /// Calculates the SequestLike Scores for all theoretical spectra.    
    let calcSequestScore scanlimits (spectrum:PeakArray<_>) scanTime chargeState isolationWindowTargetMz (theoreticalSpectra:list<TheoreticalSpectrum<Vector<float>>>) spectrumID = 
        let fCharge = float chargeState
        // measured normailzed intensity array (spectrum) minus auto-correlation
        let ms_nis =  spectrumToIntensityArrayMinusAutoCorrelation scanlimits spectrum
        // measured mass
        let ms_mass = Mass.ofMZ isolationWindowTargetMz fCharge        

        let ides = 
            theoreticalSpectra 
            |> List.fold (fun acc theoreticalSpectrum -> 
                              let lookUpResult = theoreticalSpectrum.LookUpResult
                              let sequence = lookUpResult.BioSequence
                              let p_nis = theoreticalSpectrum.TheoSpec 
                              let xcorr = calcXCorr p_nis ms_nis
                              let targetScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod true scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass sequence.Length xcorr nan nan
                              
                              let p_nis_decoy = theoreticalSpectrum.DecoyTheoSpec
                              let xcorr_decoy = calcXCorr p_nis_decoy  ms_nis
                              let decoyScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod false scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass sequence.Length xcorr_decoy nan nan
                              targetScore::decoyScore::acc  
                     ) []

        calcNormDeltaBestToRest (ides  |> List.sortBy (fun sls -> - sls.Score))        
        |> calcNormDeltaNext
        
    /// Calculates the sequestLike Scores for all theoretical spectra. Implemented using Async parallel. 
    let calcSequestScoreParallel scanlimits (spectrum:PeakArray<_>) scanTime chargeState isolationWindowTargetMz (theoreticalSpectra:list<TheoreticalSpectrum<Vector<float>>>) spectrumID = 
        let fCharge = float chargeState
        // measured normailzed intensity array (spectrum) minus auto-correlation
        let ms_nis =  spectrumToIntensityArrayMinusAutoCorrelation scanlimits spectrum
        // measured mass
        let ms_mass = Mass.ofMZ isolationWindowTargetMz fCharge     
        let computeScores theoreticalSpectrum =
            let lookUpResult = theoreticalSpectrum.LookUpResult
            let sequence = lookUpResult.BioSequence
            let p_nis = theoreticalSpectrum.TheoSpec 
            let xcorr = calcXCorr p_nis ms_nis
            let targetScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod true scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass sequence.Length xcorr nan nan
                              
            let p_nis_decoy      = theoreticalSpectrum.DecoyTheoSpec
            let xcorr_decoy      = calcXCorr p_nis_decoy ms_nis
            let decoyScore = createSearchEngineResult SearchEngineResult.SearchEngine.SEQUESTLike spectrumID lookUpResult.ModSequenceID lookUpResult.PepSequenceID lookUpResult.GlobalMod false scanTime lookUpResult.StringSequence chargeState isolationWindowTargetMz ms_mass lookUpResult.Mass sequence.Length xcorr_decoy nan nan
            [targetScore;decoyScore]
        
        let ides = 
            [for i in theoreticalSpectra ->  async {return computeScores i}]
            |> Async.Parallel 
            |> Async.RunSynchronously
            |> List.concat    
        calcNormDeltaBestToRest (ides  |> List.sortBy (fun sls -> - sls.Score))        
        |> calcNormDeltaNext

<|start_filename|>src/BioFSharp.Mz.Vis/ProteinInferenceVis.fs<|end_filename|>
namespace BioFSharp.Mz.Vis

open Plotly.NET
open BioFSharp.Mz.ProteinInference

module ProteinInference =

    let qValueHitsVisualization bandwidth (inferredProteinClassItemQValue: InferredProteinClassItemQValue[]) path (groupFiles: bool) =
        let decoy, target =
            inferredProteinClassItemQValue
            |> Array.partition (fun x -> x.InfProtClassItem.DecoyHasBetterScore)
        // Histogram with relative abundance
        let relFreqTarget = 
            FSharp.Stats.Distributions.Frequency.create bandwidth (target |> Array.map (fun x -> x.InfProtClassItem.TargetScore))
            |> Map.toArray
            |> Array.map (fun x -> fst x, (float (snd x)) / (float target.Length))
        let relFreqDecoy  =
            FSharp.Stats.Distributions.Frequency.create bandwidth (decoy |> Array.map (fun x -> x.InfProtClassItem.DecoyScore))
            |> Map.toArray
            |> Array.map (fun x -> fst x, (float (snd x)) / (float target.Length))
        // Histogram with absolute values
        let absFreqTarget =
            FSharp.Stats.Distributions.Frequency.create bandwidth (target |> Array.map (fun x -> x.InfProtClassItem.TargetScore))
            |> Map.toArray
        let absFreqDecoy  =
            FSharp.Stats.Distributions.Frequency.create bandwidth (decoy |> Array.map (fun x -> x.InfProtClassItem.DecoyScore))
            |> Map.toArray
        let histogram =
            [
                Chart.Column relFreqTarget
                |> Chart.withTraceName "Target"
                |> Chart.withAxisAnchor(Y=1);
                Chart.Column relFreqDecoy
                |> Chart.withTraceName "Decoy"
                |> Chart.withAxisAnchor(Y=1);
                Chart.Column absFreqTarget
                |> Chart.withAxisAnchor(Y=2)
                |> Chart.withMarkerStyle (Opacity = 0.)
                |> Chart.withTraceName (Showlegend = false);
                Chart.Column absFreqDecoy
                |> Chart.withAxisAnchor(Y=2)
                |> Chart.withMarkerStyle (Opacity = 0.)
                |> Chart.withTraceName (Showlegend = false)
            ]
            |> Chart.Combine

        let sortedQValues =
            inferredProteinClassItemQValue
            |> Array.map (fun x ->
                if x.InfProtClassItem.Decoy then
                    x.InfProtClassItem.DecoyScore, x.QValue
                else
                    x.InfProtClassItem.TargetScore, x.QValue
            )
            |> Array.sortBy fst

        [
            Chart.Point sortedQValues
            |> Chart.withTraceName "Q-Values";
            histogram
        ]
        |> Chart.Combine
        |> Chart.withY_AxisStyle("Relative Frequency / Q-Value",Side=StyleParam.Side.Left,Id=1, MinMax = (0., 1.))
        |> Chart.withY_AxisStyle("Absolute Frequency",Side=StyleParam.Side.Right,Id=2,Overlaying=StyleParam.AxisAnchorId.Y 1, MinMax = (0., float target.Length))
        |> Chart.withX_AxisStyle "Score"
        |> Chart.withSize (900., 900.)
        |> if groupFiles then
               Chart.SaveHtmlAs (path + @"\QValueGraph")
           else
               Chart.SaveHtmlAs (path + @"_QValueGraph")


<|start_filename|>src/BioFSharp.Mz/scripts/Testing.fsx<|end_filename|>
﻿#r "netstandard"
#r @"../../../\packages\System.Data.SQLite.Core\lib\net46\System.Data.SQLite.dll"
#r @"../../../packages\BioFSharp\lib\netstandard2.0\BioFSharp.dll"
#r @"../../../bin\BioFSharp.Mz\netstandard2.0\BioFSharp.Mz.dll"
#r @"../../../packages\FSharpAux.IO\lib\netstandard2.0\FSharpAux.IO.dll"
#r @"../../../packages\FSharpAux.IO\lib\netstandard2.0\FSharpAux.dll"
open BioFSharp
open BioFSharp.Mz

let massF:(IBioItem->float) = BioItem.initMonoisoMassWithMem

let BioSeq = "EGLDVHFVDEYEK" |> BioList.ofAminoAcidString
let BioSeqH = "EGLDVHFVDEYEK" |> BioList.ofAminoAcidString |> List.map (fun x ->AminoAcids.setModifications [ModificationInfo.Table.N15] x)

BioSeqH |> BioList.toFormula |> Formula.toString
let mono = BioSeqH |> BioList.toMonoisotopicMassWith (Formula.Table.H2O |> Formula.monoisoMass)
ModificationInfo.Table.NH3_loss |> massF

let res  = Fragmentation.Series.bOfBioList massF BioSeq
let res2 = Fragmentation.Series.bOfBioList massF BioSeqH
let light = res.[0].MainPeak.Mass
let lightLoss = res.[0].DependentPeaks.[0].Mass
light - lightLoss 

let heavy = res2.[0].MainPeak.Mass

let deph20Loss = res2.[0].DependentPeaks.[0].Mass
open Elements
System.Numerics.Complex(-0.0926829268292669, 22.0592593273255)

open System

let compareComplexNumbers (c1:System.Numerics.Complex) (c2:System.Numerics.Complex) = 
    (Math.Abs(c1.Real - c2.Real) <= System.Double.Epsilon) &&
    (Math.Abs(c1.Imaginary - c2.Imaginary) <= System.Double.Epsilon);
      
let rootO = 
    match Elements.Table.O with 
    | Tri o -> o.Root
    //| _ -> 




let rootO' =
    System.Numerics.Complex(-0.0926829268292669, 22.0592593273255),
    System.Numerics.Complex(-0.0926829268292696, -22.0592593273255)
let Zn = Multi (createMulti "Zn" (Isotopes.Table.Zn64,Isotopes.Table.Zn64.NatAbundance ) (Isotopes.Table.Zn66,Isotopes.Table.Zn66.NatAbundance) (Isotopes.Table.Zn68,Isotopes.Table.Zn68.NatAbundance) [|Isotopes.Table.Zn67;Isotopes.Table.Zn70|])
open Formula
"H2O" |> parseFormulaString |> Formula.monoisoMass
open Formula
compareComplexNumbers (fst rootO) (fst rootO)
//let iso = getSinglePhiM Table.Na 1000000. 2.
//let iso = getSinglePhiM Table.H 1000000. 2.
//let iso = getSinglePhiM Table.O 1000000. 2.
//let iso = getSinglePhiM Zn 1000. 3.
open FSharpAux
let massF : IBioItem -> float= BioItem.initMonoisoMassWithMemP  
[0.. 1000]
|> List.map (fun x -> massF AminoAcids.Ala)
[0.. 10000000]
|> PSeq.map (fun x -> massF AminoAcids.Ala)
|> Array.ofSeq
open IsotopicDistribution.MIDA
/// normalizes istopic distribution probabilities to sum up to 1.
let normalizeByProbSum minProb (pols: Polynomial list) = 
    let sum = 
        pols
        |> List.filter (fun x -> x.Probability > minProb)
        |> List.sumBy (fun x -> x.Probability)
    pols 
    |> List.map (fun x -> {x with Probability = x.Probability/sum})


"LTYYTPDYVVR"
|> BioList.ofAminoAcidString
|> BioList.toFormula
//|> fun f -> Formula.lableElement f Elements.Table.N Elements.Table.Heavy.N15
//|> (IsotopicDistribution.MIDA.ofFormula (normalizeByProbSum) 0.01 (exp(-20.)) 2)
|> (IsotopicDistribution.BRAIN.ofFormula 10 )



"LTYYTPDYVVR"
|> BioList.ofAminoAcidString
|> BioList.toFormula
|> (IsotopicDistribution.BRAIN.ofFormula 25)

BioItem.monoisoMass AminoAcids.Ala 
   
   
//[ Array.iter (fun x -> if x > 0. then do yield x ) [|-50. .. 50.|]]


let f = "N100C12O8" |> parseFormulaString 
#time


open System.Data.SQLite
let cn = SearchDB.getDBConnection @"C:\Users\david\source\repos\ProteomIQon\src\ProteomIQon\Scripts\GiadaTest.db"
let p = SearchDB.getSDBParamsBy @"C:\Users\david\source\repos\ProteomIQon\src\ProteomIQon\Scripts\GiadaTest.db"
let look = SearchDB.getThreadSafePeptideLookUpFromFileBy cn p

let psm = look 1001. 1002.

let modPep = 
    psm
    |> List.find (fun x -> x.GlobalMod = 0)

modPep.Mass
modPep.BioSequence |> List.sumBy massF |> (+) (Formula.monoisoMass Formula.Table.H2O)  

let resH = Fragmentation.Series.bOfBioList massF modPep.BioSequence 
let light' = resH.[0].MainPeak.Mass
let Loss' = resH.[0].DependentPeaks.[0].Mass

light' - Loss' 


let heavy = res2.[0].MainPeak.Mass

let deph20Loss = res2.[0].DependentPeaks.[0].Mass
heavy - deph20Loss  



//let light = res.[0].MainPeak.Mass
let heavywaterl = res2.[0].DependentPeaks.[0].Mass

heavy - light
heavywaterl - heavy


let ala = AminoAcids.Ala |> massF
let aminoLoss = 
    AminoAcids.setModification ModificationInfo.Table.NH3_loss AminoAcids.Ala
    |> massF

ala - aminoLoss 

//let alaH = 
//for i = 10 to 100000 do 
//        //AminoAcids.setModification ModificationInfo.Table.N15 AminoAcids.Ala
//        AminoAcids.Ala |> massF 
    //ModificationInfo.Table.NH3_loss 
    //|> massF

   
////let aminoLossH = 
//{ModificationInfo.Table.NH3_loss with Name= "Haralddddd"; Modify =  ModificationInfo.Table.NH3_loss.Modify >> ModificationInfo.Table.N15.Modify } 
//|> massF

//let x = {ModificationInfo.Table.NH3_loss with Name="moped"; Modify =  >>  ModificationInfo.Table.N15 .Modify } 
//x|> massF


//AminoAcids.Ala = (AminoAcids.Ala |> (AminoAcids.setModification ModificationInfo.Table.NH3_loss)  )
//(AminoAcids.Ala |> (AminoAcids.setModification ModificationInfo.Table.NH3_loss)  ) = (AminoAcids.Ala |> (AminoAcids.setModification ModificationInfo.Table.NH3_loss)  )
//(AminoAcids.Ala |> (AminoAcids.setModification ModificationInfo.Table.NH3_loss)  ) = (AminoAcids.Ala |> (AminoAcids.setModification {ModificationInfo.Table.NH3_loss with Name="moped"; Modify = ModificationInfo.Table.NH3_loss.Modify >>  ModificationInfo.Table.N15 .Modify })  )


//// Modifies existing modification (md) using (md')
//let modifyModification (md':ModificationInfo.Modification) (md:ModificationInfo.Modification) =
//    {md with Name = sprintf "%s_%s" md.Name md'.Name ; Modify = md.Modify >>  md'.Modify} 



//let heavyNh3 = modifyModification ModificationInfo.Table.N15 ModificationInfo.Table.NH3_loss

//ModificationInfo.Table.NH3_loss |> massF

//heavyNh3 |> massF

//ModificationInfo.Table.N15 |> massF






<|start_filename|>src/BioFSharp.Mz/PeakList.fs<|end_filename|>
﻿namespace BioFSharp.Mz

type PeakList<[<EqualityConditionalOn; ComparisonConditionalOn >]'a when 'a :> IPeak> = list<'a>

module PeakList =
    
    /// Builds a new PeakArray whose elements are the results of applying the given function to each 
    /// of the elements of the PeakArray.
    let map f pkl : PeakList<_> = 
        List.map f pkl

    /// Iterates the mz and intensity lists and creates a Peak(mz,intensity) for each value pair. 
    /// Returns a new PeakList. 
    let zip (mz:list<float>) (intensity:list<float>) : PeakList<_> = 
        List.map2 (fun m i -> Peak(m,i)) mz intensity

    /// Iterates the PeakList and unzips the fields of each peak into two seperate lists, the first containing 
    /// the mz values, the second the intensities. 
    let unzip (pkl : PeakList<_>) = 
        let rec loop (l:PeakList<_>) mz intens =
            match l with
            | h::t -> loop t (h.Mz::mz) (h.Intensity::intens)
            | _ -> List.rev mz ,List.rev intens           // TODO: Continuation passing
        loop pkl [] []

    /// Creates a Peak(mz,intensity) for each value pair in mzIntensity. 
    /// Returns a new PeakList.
    let zipMzInt (mzIntensity:list<float*float>) : PeakList<_> = 
        List.map (fun (m,i) -> Peak(m,i)) mzIntensity


<|start_filename|>src/BioFSharp.Mz/FSharp.Stats.fs<|end_filename|>
﻿namespace BioFSharp.Mz

open FSharp.Stats
open FSharp.Stats.Signal
open System


module NonLinearRegression' = 
    ///
    let standardErrorOfPrediction dOF (predicted:float []) (actual:float [])  =
        let n = actual.Length-1 |> float 
        match n with
        | x when x > dOF -> 
            let sumOfResSq = Array.fold2 (fun acc yReal yPred  -> acc + ((yPred-yReal)**2.) ) 0.0  actual predicted
            sqrt( (sumOfResSq / (n - dOF)))
        | _             -> nan

module PeakDetection = 

    module Wavelet = 

        type WaveletPeak = {
            XAxisIndex  : int 
            ScaleIndex  : int
            Scale       : float
            Correlation : float
            }

        type Gaussian = {
            Amplitude   : float
            XLoc        : float
            Stdev       : float
            Trace       : (float*float) []
            Function    : float -> float
            }

        type PeakGroup = {   
            Start       : float
            End         : float
            Data        : (float*float)[]
            Fits        : Gaussian list
            SumTrace    : float []
            Convolved   : bool
            }

        type Parameters = {
            Borderpadding           : int 
            BorderPadMethod         : Padding.BorderPaddingMethod 
            InternalPaddingMethod   : Padding.InternalPaddingMethod 
            HugeGapPaddingMethod    : Padding.HugeGapPaddingMethod 
            HugeGapPaddingDistance  : float
            MaxPeakLength           : float 
            NoiseQuantile           : float 
            MinSNR                  : float
            }

        ///
        let private getColLowBound  (xData: float[]) centerIdx mzTol = 
            let rec loop  (xData: float []) center i mzTol =
                if i <= 0 then 0
                else                                              
                    match (xData.[center] - xData.[i]) with 
                    | x when x >= mzTol -> i+1   
                    | x when x <= mzTol -> loop xData center (i-1) mzTol
            loop xData centerIdx (centerIdx-1) mzTol

        /// 
        let private getColHighBound (xData: float []) centerIdx mzTol = 
            let rec loop (xData: float []) center i mzTol = 
                if i >= xData.Length-1 then xData.Length-1
                else                                             
                    match xData.[i] - xData.[center] with 
                    | x when x >= mzTol -> i-1
                    | x when x <= mzTol -> loop xData center (i+1) mzTol
            loop xData centerIdx (centerIdx+1) mzTol

        ///
        let private getScales maxScale averageDistance =
            let tmp =
                [|averageDistance*2. .. round 2 (averageDistance/2.) .. maxScale|]
            tmp

        ///
        let private getMinPeakDistance averageDistance = averageDistance*3.

        ///
        let private aic n nParams likelihood = 
            (float n * Math.Log(likelihood / float n)) + 2. * (float nParams) 

        ///
        let private optimizeWaveletFits (origData:(float*float) []) (groupedPeaks:(WaveletPeak*Gaussian)list) = 
            match groupedPeaks with 
            | (wp,g)::[] ->
                let min = g.XLoc - 2. * g.Stdev
                let max = g.XLoc + 2. * g.Stdev
                let data = origData |> Array.filter (fun (x,y) -> x >= min  && x <= max )
                let yHat = origData |> Array.map (fst >> g.Function)
                {Start=min;End=max;Data=data;Fits=[g];SumTrace=yHat;Convolved=false}
            | _ -> 
                let min  = groupedPeaks |> List.map (fun (wp,g) -> g.XLoc - 2. * g.Stdev) |> List.min
                let max  = groupedPeaks |> List.map (fun (wp,g) -> g.XLoc + 2. * g.Stdev) |> List.max
                let xy = origData |> Array.filter (fun (x,y) -> x >= min && x <= max) 
                let x,y = xy |> Array.unzip        
                let bestGroup = 
                    groupedPeaks
                    |> List.map (fun (wp,g) -> 
                        let yHat = x |> Array.map g.Function
                        g,yHat
                    )
                    |> FSharpAux.List.powerSetOf
                    |> List.minBy (fun groupedPeaks ->    
                        let m = Matrix.ofJaggedArray (groupedPeaks |> List.map snd |> Array.ofList)
                        let yHat'= m |> Matrix.mapiCols (fun i v -> v |> Vector.sum)
                        let rss = Seq.map2 (fun x y -> (x-y)**2. ) y yHat' |> Seq.sum
                        aic m.NumCols m.NumRows rss
                    )
                    |> List.map fst
                let sumTrace = 
                    let m = Matrix.ofJaggedArray (bestGroup |> List.map (fun g ->  origData |> Array.map (fst >> g.Function)) |> Array.ofList )
                    let yHat' = m |> Matrix.mapiCols (fun i v -> v |> Vector.sum)
                    yHat'
                    |> Array.ofSeq            
                {Start=min;End=max;Data=xy;Fits=bestGroup;SumTrace=sumTrace;Convolved=true}
             
        ///        
        let identifyPeaksBy (borderpadding:int) (borderPadMethod:Padding.BorderPaddingMethod) (internalPaddingMethod:Padding.InternalPaddingMethod) (hugeGapPaddingMethod:Padding.HugeGapPaddingMethod) 
            (maxDistance:float) maxPeakLength noiseQuantile minSNR (trace:(float*float)[]) =
            let maxScale = maxPeakLength / 6.
            let averageDistance = Padding.HelperFunctions.getMedianSpacing trace (-)
            let minDistancePadding = averageDistance / 2. 
            let paddedData = Padding.pad trace minDistancePadding maxDistance (-) (+) borderpadding borderPadMethod internalPaddingMethod hugeGapPaddingMethod
            let minPeakDistance = getMinPeakDistance minDistancePadding
            let scales = getScales maxScale minDistancePadding
            let transformations = 
                scales
                |> Array.map Wavelet.createRicker
                |> Array.map (fun x -> ContinuousWavelet.transform paddedData (-) borderpadding x |> Array.unzip)
            let xVals = fst transformations.[0]
            let corrMatrix = 
                transformations
                |> Array.map snd
                |> Matrix.ofJaggedArray
            let noiseLevel = (corrMatrix.Row 0).ToArray() |> FSharp.Stats.Quantile.compute noiseQuantile
            let peakMatrix = 
                corrMatrix
                |> Matrix.mapi (fun m n x -> 
                    if n < 2 || n >= corrMatrix.NumCols-3 then 0. 
                    elif x > corrMatrix.[m,n-1] && x > corrMatrix.[m,n-2] && x > corrMatrix.[m,n+1] && x > corrMatrix.[m,n+2] && x >= (minSNR*noiseLevel) then 
                        x
                    else 0.
                )
            let maxCorrScale = 
                peakMatrix 
                |> Matrix.mapiCols (fun i x ->
                    let mutable maxScaleIdx,maxScale,maxCorr = 0,0.,0.
                    for scale = 1 to x.Length-1 do 
                        if x.[scale] > maxCorr then 
                            maxScaleIdx <- scale 
                            maxScale    <- scales.[scale]
                            maxCorr     <- x.[scale]
                    {XAxisIndex=i ;ScaleIndex=maxScaleIdx; Scale=maxScale; Correlation=maxCorr}   
                )
                |> Array.ofSeq
            let mergedPeaks = 
                let finalPeaks = ResizeArray<WaveletPeak>(100)
                for i = 2 to maxCorrScale.Length-3 do    
                    let currentPeak = maxCorrScale.[i]
                    let currentCorr = currentPeak.Correlation
                    if currentCorr < maxCorrScale.[i-1].Correlation ||
                        currentCorr < maxCorrScale.[i-2].Correlation ||
                        currentCorr < maxCorrScale.[i+1].Correlation ||
                        currentCorr < maxCorrScale.[i+2].Correlation ||
                        currentCorr = 0. then ()
                    else 
                        if currentPeak.ScaleIndex > 0 (*nouse*) then 
                            let lowBound  = getColLowBound xVals i (minPeakDistance/2.)
                            let highBound = getColHighBound xVals i (minPeakDistance/2.)
                            let mutable maxCorr = 0.0
                            let mutable maxCol = 0
                            for j = lowBound to highBound do 
                                let scale = maxCorrScale.[j].ScaleIndex
                                if corrMatrix.[scale,j] > maxCorr && scale > 0  then
                                    maxCorr <- corrMatrix.[scale,j]
                                    maxCol <- j 
                            let refinedPeak = maxCorrScale.[maxCol] 
                            if finalPeaks.Count = 0 then 
                                finalPeaks.Add refinedPeak
                            else 
                                let prevPeak = finalPeaks.[finalPeaks.Count-1]
                                let mzDiff = xVals.[refinedPeak.XAxisIndex] - xVals.[prevPeak.XAxisIndex]
                                if mzDiff > minPeakDistance then
                                    finalPeaks.Add refinedPeak
                                elif refinedPeak.Correlation > prevPeak.Correlation then
                                    finalPeaks.[finalPeaks.Count-1] <- refinedPeak
                finalPeaks.ToArray()
            let refinedPeaks = 
                mergedPeaks 
                |> Array.map (fun wp -> 
                    let loc,apex = paddedData |> Array.minBy (fun x -> abs(fst x - (xVals.[wp.XAxisIndex])))
                    let f = FSharp.Stats.Fitting.NonLinearRegression.Table.gaussModel.GetFunctionValue ([|apex;loc;wp.Scale|] |> vector)
                    let peakSim = 
                        trace 
                        |> Array.map fst 
                        |> Array.map (fun x -> x,f x)
                    let gaussian = 
                        {
                            Amplitude   = apex
                            XLoc        = loc
                            Stdev       = wp.Scale
                            Trace       = peakSim
                            Function = f
                        }
                    wp,gaussian 
                )
                |> Array.sortBy (fun (wp,g) -> wp.XAxisIndex)
                |> List.ofArray
                |> List.filter (fun (wp,g) -> g.Stdev > 0. && g.Amplitude > 0.)
            let groupedPeaks = 
                let isOverlapping (currentPeak:(WaveletPeak*Gaussian)) (candidate:(WaveletPeak*Gaussian)) = 
                    (((snd candidate).XLoc < ((snd currentPeak).XLoc + ((snd currentPeak).Stdev * 2.))) && ((snd candidate).XLoc > ((snd currentPeak).XLoc - ((snd currentPeak).Stdev * 2.))) ) ||
                    ( ((snd currentPeak).XLoc < ((snd candidate).XLoc + ((snd candidate).Stdev * 2.))) && ((snd currentPeak).XLoc > ((snd candidate).XLoc - ((snd candidate).Stdev * 2.)))     ) 
                let overlappingPeaks currentFamily peaksLeft = 
                    let rec findOverlappingPeaks  (overlappingPeaks:(WaveletPeak*Gaussian) list) (nonOverlappingPeaks:(WaveletPeak*Gaussian) list) (currentPeak:(WaveletPeak*Gaussian)) (peaksLeft:(WaveletPeak*Gaussian) list) =
                        match peaksLeft with 
                        | [] -> overlappingPeaks,nonOverlappingPeaks
                        | h::t ->
                            if isOverlapping currentPeak h then 
                                findOverlappingPeaks (h::overlappingPeaks) nonOverlappingPeaks currentPeak t 
                            else 
                                findOverlappingPeaks overlappingPeaks (h::nonOverlappingPeaks) currentPeak t 
                              
                    let rec loop newPeakFamily currentFamily peaksLeft =
                        match currentFamily with 
                        | [] -> newPeakFamily,peaksLeft
                        | cp::t -> 
                            let overlappingPeaks,nonOverlappingPeaks = findOverlappingPeaks [] [] cp peaksLeft
                            loop (overlappingPeaks@newPeakFamily) t nonOverlappingPeaks
                    loop [] currentFamily peaksLeft             
                let acc = ResizeArray<(WaveletPeak*Gaussian)list>()
                let rec loop (acc:ResizeArray<(WaveletPeak*Gaussian)list>) currentFam peaksLeft =
                    match peaksLeft with 
                    | [] -> 
                        currentFam 
                        |> acc.Add 
                        acc |> Seq.toList
                    | cP::t ->
                        match currentFam with 
                        | [] -> 
                            loop acc (cP::currentFam) (t)                     
                        | cf  ->
                            let overlapping = overlappingPeaks cf (cP::t) 
                            match overlapping with 
                            | [],rest   -> 
                                cf 
                                |> acc.Add 
                                loop acc [] rest
                            | newPeaks,rest -> 
                                let cf' = newPeaks@cf
                                loop acc cf' rest 
                loop acc [] refinedPeaks                     
            let origData = 
                xVals 
                |> Array.map (fun x -> 
                    paddedData |> Array.minBy (fun (xx,yy) -> abs(xx-x) )
                )

            let fits = 
                groupedPeaks
                |> List.filter  (List.isEmpty >> not)
                |> List.map (optimizeWaveletFits origData)
            let final =
                let m = 
                    fits 
                    |> List.map (fun x -> x.SumTrace) 
                    |> Array.ofList
                m 
                |> JaggedArray.transpose
                |> Array.map (fun v -> 
                    let max = v |> Seq.max
                    v 
                    |> Array.map (fun x -> if x < max then 0. else x)
                )
                |> JaggedArray.transpose
                |> Array.mapi (fun i trace -> 
                    let g = fits.[i]
                    let y' = trace
                    let max = y' |> Array.max
                    let sumTrace' = 
                        Array.map2 (fun (x,y) y' -> x,y,y') origData y'
                        |> Array.choose (fun (x,y,y') -> if y' > 0. && y' > max * 0.01 then Some (x,y) else None)           
                    if Array.isEmpty sumTrace' then 
                        None
                    else
                        let start' =  sumTrace' |> Array.minBy fst |> fst
                        let end'   =  sumTrace' |> Array.maxBy fst |> fst
                        Some {g with Start=start';End=end';Data=sumTrace' }
                )
                |> Array.choose id
                |> Array.map (fun pk -> 
                                let x,y = pk.Data |> Array.unzip
                                let traces = 
                                    Array.map2 (fun x y -> 
                                        let yHats = 
                                            pk.Fits 
                                            |> List.mapi (fun i g -> i,g.Function x)
                                        let maxIndex,value = List.maxBy snd yHats
                                        yHats 
                                        |> List.map (fun (i,value) -> if i = maxIndex then x,y else x,0.)
                                        |> Array.ofList
                                    ) x y
                                    |> JaggedArray.transpose
                            
                                let refined = 
                                    pk.Fits
                                    |> List.mapi (fun i g -> 
                                        let trace = 
                                            let tmp = 
                                                traces.[i] |> Array.filter (fun (x,y) -> y > 0.)
                                            if Array.isEmpty tmp then 
                                                None
                                            else    
                                                let maxX,maxY = tmp |> Array.maxBy snd 
                                                tmp 
                                                |> Array.filter (fun x -> snd x > maxY * 0.01)
                                                |> Some 
                                        if trace.IsNone then None 
                                        else
                                        Some {g with Trace = trace.Value}
                                    )
                                    |> List.choose id
                                {pk with Fits = refined}
                            )
                |> List.ofSeq
            final
           
        ///        
        let identifyPeaks (parameters:Parameters) (trace:(float*float)[]) =
            identifyPeaksBy parameters.Borderpadding parameters.BorderPadMethod parameters.InternalPaddingMethod parameters.HugeGapPaddingMethod
                parameters.HugeGapPaddingDistance parameters.MaxPeakLength parameters.NoiseQuantile parameters.MinSNR trace

        //
        let substractBaseLine (yData:float []) =
            let baseLine = FSharp.Stats.Signal.Baseline.baselineAls 10 6 0.05 yData |> Array.ofSeq
            Array.map2 (fun y b ->
                           let c = y - b
                           if c < 0. then 0. else c
                       ) yData baseLine


<|start_filename|>docs/tools/BioFSharp.Mz.fsx<|end_filename|>
(**
---
title: BioFSharp.Mz
category: Description
categoryindex: 1
index: 1
---
*)

<|start_filename|>src/BioFSharp.Mz/Quantification.fs<|end_filename|>
﻿namespace BioFSharp.Mz

module Quantification =

    open System
    open SignalDetection
    open FSharp.Stats
    open Fitting.NonLinearRegression
    open Signal.PeakDetection


    module Integration = 

        /// Returns the estimated area beneath the data using the trapezoidal rule. Requires uniform spaced datapoints. 
        let trapezEstAreaOfUniform (xData: float[]) (yData: float[]) =
            if xData.Length <> yData.Length then
                failwith "Both input collections must be of the same size"
            // lower border of area estimation
            let a = xData.[0]
            // upper border of area estimation
            let b = xData.[xData.Length-1]
            let mutable trapezSum = 0.0
            for i = 1 to xData.Length-2 do
                trapezSum <- trapezSum + 2. * yData.[i]
            trapezSum <- trapezSum + yData.[0]
            trapezSum <- trapezSum + yData.[xData.Length-1]
            let area = ( (b-a) / (2. * float xData.Length) ) * trapezSum
            area    
         
        /// Returns the estimated area beneath the data using the trapezoidal rule.
        let trapezEstAreaOf (xData: float[]) (yData: float[]) =
            if xData.Length <> yData.Length then
                failwith "Both input collections must be of the same size"
            let mutable trapezSum = 0.0
            for i = 0 to xData.Length-2 do
                trapezSum <- trapezSum + (xData.[i+1] - xData.[i]) * (yData.[i+1] + yData.[i])
            let area = 0.5 * trapezSum
            area    

    module ParameterEstimation =

        ///
        type CaruanaEstimates = {
            Amplitude: float
            MeanX    : float
            STD      : float
            FWHM     : float
            }

        /// 
        let createCaruanaEstimates amplitude meanX std fwhm = {
            Amplitude=amplitude; MeanX=meanX; STD=std; FWHM=fwhm }

        /// Returns the value of the standard deviation
        let toSTD fwhm = fwhm / ( 2. * sqrt(2. * log 2.) )

        /// Returns the full width at half maximum  
        let toFWHM std = 2. * sqrt(2. * log 2.) * std

        /// Returns the yValue of a gauss function at a given position x.
        let gaussFunc amplitude meanX std x = 
            amplitude * exp((-1.)*((((x-meanX)**2.)/(2.*std**2.))))
  
        /// Estimates the Parameters of a Gaussian function
        /// Warning: This method is sensitive to noisy data. If the noise level of the input parameters is high, smoothing of 
        /// the data is strongly recommended. 
        let caruanaAlgorithm (mzData:float []) (intensityData:float []) =
            let mzData,intensityData = 
                Array.zip mzData intensityData
                |> Array.filter (fun (_,intensity) -> intensity <> 0.)
                |> Array.unzip
                |> fun (x,y) -> Vector.ofArray x, y
            if mzData.Length < 3 || intensityData.Length < 3 then Option.None 
            else 
            let logTransIntensityData = 
                intensityData
                |> Array.map (fun x -> 
                                if x <= 0. then log 1. 
                                else log x
                             )
                |> Vector.ofArray                          
            let polCoeff = FSharp.Stats.Fitting.LinearRegression.OrdinaryLeastSquares.Polynomial.coefficient 2 mzData logTransIntensityData 
            // f(x) = a1 + a2 * x + a3 * x**2
            let a = polCoeff.[0]
            let b = polCoeff.[1]
            let c = polCoeff.[2]
            let amplitude = exp(a-((b**2.)/(4.*c)))
            let meanX = -b/(2.*c) 
            let fwhm  = sqrt(-1./(2.*c)) * sqrt(2.*log(2.))*2.
            let std   = toSTD fwhm
            Some (createCaruanaEstimates amplitude meanX std fwhm)


        ///
        type EstimatedMoments = {
            ModeY : float 
            MeanX : float 
            Var   : float 
            Std   : float 
            Skew  : float 
            Tau   : float
            }

        ///
        let createEstimatedMoments modeY meanX var std skew tau = {
            ModeY = modeY 
            MeanX = meanX 
            Var   = var   
            Std   = std   
            Skew  = skew  
            Tau   = tau      
            }

        ///
        let meanOfGaussian xData yData = 
            Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y*x) ) (0.,0.)  xData yData 
            |> fun (sum,acc) -> acc / sum

        ///
        let varianceBy mean xData yData =     
            Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y * (x - mean)**2.)  ) (0.,0.)  xData yData 
            |> fun (sum,acc) -> acc/sum

        ///
        let varianceOf xData yData =     
            let mean = meanOfGaussian xData yData
            varianceBy mean xData yData 

        ///
        let skewBy mean var xData yData = 
            let m3   =  
                Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y * (x - mean)**3.) ) (0.,0.)  xData yData 
                |> fun (sum,acc) -> acc/sum
            m3 / (var**1.5)

        ///
        let skewOf xData yData = 
            let mean = meanOfGaussian xData yData
            let var  = varianceBy mean xData yData
            skewBy mean var xData yData 
    
        ///    
        let estTau stdev skew  =
            stdev * (skew/2.) ** 0.333

        ///
        let methodOfMoments (p: IdentifiedPeak) =
            try
            let modeY   = p.Apex.YVal
            let meanX   = meanOfGaussian p.XData p.YData
            let var     = varianceOf p.XData p.YData
            let std     = sqrt(var)
            let skew    = skewBy meanX var p.XData p.YData
            let tau     = estTau std skew
            Some (createEstimatedMoments modeY meanX var std skew tau)
            with 
            | _ -> Option.None
        ///
        let estimatePeakIntegrity (p: IdentifiedPeak) = 
            let zipped = Array.zip p.XData p.YData 
            let left  = zipped |> Array.exists (fun (x,y) -> x < p.Apex.XVal && y > 0.2 * p.Apex.YVal) 
            let right = zipped |> Array.exists (fun (x,y) -> x > p.Apex.XVal && y > 0.2 * p.Apex.YVal)
            left && right
            
        ///
        let estimateMoments (p: IdentifiedPeak ) =
            if estimatePeakIntegrity p then 
                methodOfMoments p 
            else 
                let xData,yData = 
                    Array.zip p.XData p.YData 
                    |> Array.filter (fun (x,y) -> y > 0.3 * p.Apex.YVal)
                    |> Array.unzip
                match caruanaAlgorithm xData yData with 
                | Some caruanaEst -> 
                    let skew = skewBy caruanaEst.MeanX  (caruanaEst.STD**2.) p.XData p.YData
                    let tau  = estTau caruanaEst.STD skew
                    Some (createEstimatedMoments p.Apex.YVal caruanaEst.MeanX (caruanaEst.STD**2.) caruanaEst.STD skew tau)
                | Option.None            -> 
                    methodOfMoments p 
      
 
    module HULQ = 

        open Integration
        open ParameterEstimation
        open LevenbergMarquardtConstrained

        type PeakModel = 
            | Gaussian of Model
            | EMG of Model 

        type FittedPeak = {
            Model                       : PeakModel 
            EstimatedParams             : float []
            StandardErrorOfPrediction   : float
            YPredicted                  : float []
            Area                        : float
            }

        ///
        let createFittedPeak model estimatedParams standardErrorOfPrediction yPredicted area = {
            Model                       = model                       
            EstimatedParams             = estimatedParams            
            StandardErrorOfPrediction   = standardErrorOfPrediction  
            YPredicted                  = yPredicted                 
            Area                        = area
            }

        type QuantifiedPeak = {
            Model                       : PeakModel option 
            YPredicted                  : float []
            EstimatedParams             : float []
            StandardErrorOfPrediction   : float
            Area                        : float
            MeasuredApexIntensity       : float
            }

        ///
        let createQuantifiedPeak model yPredicted estimatedParams standardErrorOfPrediction area measuredApexIntensity = {
            Model                       = model                       
            YPredicted                  = yPredicted                  
            EstimatedParams             = estimatedParams             
            StandardErrorOfPrediction   = standardErrorOfPrediction   
            Area                        = area      
            MeasuredApexIntensity       = measuredApexIntensity
            }

        /// Return Option
        let getPeakBy (peaks:IdentifiedPeak []) x =
            let isPartOfPeak x (p:IdentifiedPeak) = p.LeftEnd.XVal < x && p.RightEnd.XVal > x 
            match Array.tryFind (isPartOfPeak x) peaks with
            | Some p        ->  p
            | Option.None   -> Array.minBy (fun p -> abs(p.Apex.XVal - x)) peaks 

        ///
        let integralOfGaussian (coeffs:vector ) =
            sqrt(2.) * coeffs.[0] * coeffs.[2] * sqrt(Math.PI)

        ///
        let integralOfEMGBy meanX sigma tau intensityAtx xPos =    
            let fst = intensityAtx * (sqrt(Math.PI)*tau)
            let intg =         
                let firstT = sqrt(Math.PI) / 2.
                let errorF = 
                    ((xPos-meanX)/(sqrt(2.)*sigma)) - (sigma/(sqrt(2.)*tau))
                    |> FSharp.Stats.SpecialFunctions.Errorfunction.Erf
                firstT * (1.+errorF)
            let exp = exp ( 0.5*((sigma/tau)**2.) - ((xPos-meanX)/tau) )
            fst / (intg * exp)

        ///
        let calcArea model (estParams:vector) =
            match model with 
            | Gaussian f -> 
                integralOfGaussian estParams
            | EMG      f -> 
                let intensityAtx = f.GetFunctionValue estParams estParams.[1]
                integralOfEMGBy estParams.[1] estParams.[2] estParams.[3] intensityAtx estParams.[1]

        ///
        let tryFit model solverOptions lowerBound upperBound xData yData =
            let modelF = 
                match model with 
                | Gaussian f -> f
                | EMG f      -> f       
            try
                let estParams = estimatedParamsVerbose modelF solverOptions 0.001 10.0 lowerBound upperBound xData yData |> Seq.last                                  
                let y' = Array.map (fun xValue -> modelF.GetFunctionValue estParams xValue) xData
                let sEoE_Gauss = NonLinearRegression'.standardErrorOfPrediction (float solverOptions.InitialParamGuess.Length) y' yData
                let area = calcArea model estParams
                if nan.Equals area || infinity.Equals area || infinity.Equals (area * -1.) || nan.Equals sEoE_Gauss then 
                    Option.None 
                else
                    Option.Some (createFittedPeak model (estParams.ToArray()) sEoE_Gauss y' area)
            with 
            | ex -> Option.None      

        ///
        let tryFitGaussian initAmp initMeanX initStdev xData yData =
            let deltaAmp = initAmp * 0.3
            let deltaX = Math.Min((Array.max xData - Array.min xData) / 2.,10.) 
            let deltaStabw = initStdev * 0.5
            let lowerBound = [|Math.Max(initAmp-deltaAmp,0.); Math.Max(initMeanX-deltaX,0.); Math.Max(initStdev-deltaStabw,0.)|] |> vector
            let upperBound = [|initAmp+deltaAmp; initMeanX+deltaX; initStdev+deltaStabw|] |> vector
            let solverOptions = createSolverOption 0.00001 0.001 1000 [|initAmp;initMeanX;initStdev|]
            tryFit (Gaussian Table.gaussModel) solverOptions lowerBound upperBound xData yData

        ///
        let tryFitEMG initAmp initMeanX initStdev initTau xData yData =
            let deltaAmp = initAmp * 0.3
            let deltaX = Math.Min((Array.max xData - Array.min xData) / 2.,10.) 
            let deltaStabw = initStdev * 0.5
            let deltaTau = initTau
            let lowerBound = [|Math.Max(initAmp-deltaAmp,0.); Math.Max(initMeanX-deltaX,0.); Math.Max(initStdev-deltaStabw,0.);Math.Max(initTau-deltaTau,0.00001)|] |> vector
            let upperBound = [|initAmp+deltaAmp; initMeanX+deltaX; initStdev+deltaStabw; initTau+deltaTau|] |> vector
            let solverOptions = createSolverOption 0.00001 0.001 1000 [|initAmp;initMeanX;initStdev;initTau|]
            tryFit (EMG Table.emgModel) solverOptions lowerBound upperBound xData yData

        ///
        let selectModel (models: FittedPeak []) = 
            Array.minBy (fun (x:FittedPeak) -> x.StandardErrorOfPrediction) models


        ///
        let quantifyPeak (p: IdentifiedPeak ) =
            match estimateMoments p with 
            | Some moments -> 
                let gaussian = tryFitGaussian moments.ModeY moments.MeanX moments.Std p.XData p.YData
                let emg      = tryFitEMG moments.ModeY moments.MeanX moments.Std moments.Tau p.XData p.YData
                match gaussian, emg with 
                | Some g, Some emg -> 
                    let finalFit   = selectModel [|g;emg|]
                    createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction finalFit.Area p.Apex.YVal
                | Some finalFit, Option.None  | Option.None, Some finalFit -> 
                    createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction finalFit.Area  p.Apex.YVal
                | _                       ->
                    let area = trapezEstAreaOf p.XData p.YData
                    createQuantifiedPeak Option.None [||] [||] nan area p.Apex.YVal
            | Option.None   -> 
                let area = trapezEstAreaOf p.XData p.YData
                createQuantifiedPeak Option.None [||] [||] nan area p.Apex.YVal              

          

<|start_filename|>tests/BioFSharp.Mz.Tests/Main.fs<|end_filename|>
﻿module BioFSharp.Mz.Tests
open Expecto

[<EntryPoint>]
let main argv =
    Tests.runTestsInAssembly defaultConfig argv


<|start_filename|>src/BioFSharp.Mz/Quantification - Copy.fs<|end_filename|>
﻿namespace BioFSharp.Mz

//module Quantification =

//    open System
//    open SignalDetection
//    open FSharp.Stats
//    open Fitting.NonLinearRegression
//    open Signal.PeakDetection
//    open PeakDetection.Wavelet

//    type PeakToQuantify<'T> = {
//        ApexY : float 
//        ApexX : float
//        XData : float []
//        YData : float []
//        ParameterEstimates: 'T option
//        }

//    type SignalEntity<'T> = {
//        Start                   : float 
//        End                     : float
//        XData                   : float []
//        YData                   : float []
//        PeaksToQuantify         : seq<PeakToQuantify<'T>>
//        HasParameterEstimates   : bool
//        HasMultiplePeaks        : bool
//    }

//    module Integration = 

//        /// Returns the estimated area beneath the data using the trapezoidal rule. Requires uniform spaced datapoints. 
//        let trapezEstAreaOfUniform (xData: float[]) (yData: float[]) =
//            if xData.Length <> yData.Length then
//                failwith "Both input collections must be of the same size"
//            // lower border of area estimation
//            let a = xData.[0]
//            // upper border of area estimation
//            let b = xData.[xData.Length-1]
//            let mutable trapezSum = 0.0
//            for i = 1 to xData.Length-2 do
//                trapezSum <- trapezSum + 2. * yData.[i]
//            trapezSum <- trapezSum + yData.[0]
//            trapezSum <- trapezSum + yData.[xData.Length-1]
//            let area = ( (b-a) / (2. * float xData.Length) ) * trapezSum
//            area    
         
//        /// Returns the estimated area beneath the data using the trapezoidal rule.
//        let trapezEstAreaOf (xData: float[]) (yData: float[]) =
//            if xData.Length <> yData.Length then
//                failwith "Both input collections must be of the same size"
//            let mutable trapezSum = 0.0
//            for i = 0 to xData.Length-2 do
//                trapezSum <- trapezSum + (xData.[i+1] - xData.[i]) * (yData.[i+1] + yData.[i])
//            let area = 0.5 * trapezSum
//            area    

//    module ParameterEstimation =

//        ///
//        type CaruanaEstimates = {
//            Amplitude: float
//            MeanX    : float
//            STD      : float
//            FWHM     : float
//            }

//        /// 
//        let createCaruanaEstimates amplitude meanX std fwhm = {
//            Amplitude=amplitude; MeanX=meanX; STD=std; FWHM=fwhm }

//        /// Returns the value of the standard deviation
//        let toSTD fwhm = fwhm / ( 2. * sqrt(2. * log 2.) )

//        /// Returns the full width at half maximum  
//        let toFWHM std = 2. * sqrt(2. * log 2.) * std

//        /// Returns the yValue of a gauss function at a given position x.
//        let gaussFunc amplitude meanX std x = 
//            amplitude * exp((-1.)*((((x-meanX)**2.)/(2.*std**2.))))
  
//        /// Estimates the Parameters of a Gaussian function
//        /// Warning: This method is sensitive to noisy data. If the noise level of the input parameters is high, smoothing of 
//        /// the data is strongly recommended. 
//        let caruanaAlgorithm (mzData:float []) (intensityData:float []) =
//            let mzData,intensityData = 
//                Array.zip mzData intensityData
//                |> Array.filter (fun (_,intensity) -> intensity <> 0.)
//                |> Array.unzip
//                |> fun (x,y) -> Vector.ofArray x, y
//            if mzData.Length < 3 || intensityData.Length < 3 then Option.None 
//            else 
//            let logTransIntensityData = 
//                intensityData
//                |> Array.map (fun x -> 
//                                if x <= 0. then log 1. 
//                                else log x
//                             )
//                |> Vector.ofArray                          
//            let polCoeff = FSharp.Stats.Fitting.LinearRegression.OrdinaryLeastSquares.Polynomial.coefficient 2 mzData logTransIntensityData 
//            // f(x) = a1 + a2 * x + a3 * x**2
//            let a = polCoeff.[0]
//            let b = polCoeff.[1]
//            let c = polCoeff.[2]
//            let amplitude = exp(a-((b**2.)/(4.*c)))
//            let meanX = -b/(2.*c) 
//            let fwhm  = sqrt(-1./(2.*c)) * sqrt(2.*log(2.))*2.
//            let std   = toSTD fwhm
//            Some (createCaruanaEstimates amplitude meanX std fwhm)


//        ///
//        type EstimatedMoments = {
//            ModeY : float 
//            MeanX : float 
//            Var   : float 
//            Std   : float 
//            Skew  : float 
//            Tau   : float
//            }

//        ///
//        let createEstimatedMoments modeY meanX var std skew tau = {
//            ModeY = modeY 
//            MeanX = meanX 
//            Var   = var   
//            Std   = std   
//            Skew  = skew  
//            Tau   = tau      
//            }

//        ///
//        let meanOfGaussian xData yData = 
//            Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y*x) ) (0.,0.)  xData yData 
//            |> fun (sum,acc) -> acc / sum

//        ///
//        let varianceBy mean xData yData =     
//            Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y * (x - mean)**2.)  ) (0.,0.)  xData yData 
//            |> fun (sum,acc) -> acc/sum

//        ///
//        let varianceOf xData yData =     
//            let mean = meanOfGaussian xData yData
//            varianceBy mean xData yData 

//        ///
//        let skewBy mean var xData yData = 
//            let m3   =  
//                Array.fold2 (fun (sum,acc) x y -> sum+y, acc + (y * (x - mean)**3.) ) (0.,0.)  xData yData 
//                |> fun (sum,acc) -> acc/sum
//            m3 / (var**1.5)

//        ///
//        let skewOf xData yData = 
//            let mean = meanOfGaussian xData yData
//            let var  = varianceBy mean xData yData
//            skewBy mean var xData yData 
    
//        ///    
//        let estTau stdev skew  =
//            stdev * (skew/2.) ** 0.333

//        ///
//        let methodOfMoments (p: PeakToQuantify<'a>) =
//            try
//            let modeY   = p.ApexY
//            let meanX   = meanOfGaussian p.XData p.YData
//            let var     = varianceOf p.XData p.YData
//            let std     = sqrt(var)
//            let skew    = skewBy meanX var p.XData p.YData
//            let tau     = estTau std skew
//            Some (createEstimatedMoments modeY meanX var std skew tau)
//            with 
//            | _ -> Option.None
//        ///
//        let estimatePeakIntegrity (p: PeakToQuantify<'a>) = 
//            let zipped = Array.zip p.XData p.YData 
//            let left  = zipped |> Array.exists (fun (x,y) -> x < p.ApexX && y > 0.2 * p.ApexY) 
//            let right = zipped |> Array.exists (fun (x,y) -> x > p.ApexX && y > 0.2 * p.ApexY)
//            left && right

//        ///
//        let estimateMoments (p: PeakToQuantify<'a> ) =
//            if estimatePeakIntegrity p then 
//                methodOfMoments p 
//            else 
//                let xData,yData = 
//                    Array.zip p.XData p.YData 
//                    |> Array.filter (fun (x,y) -> y > 0.3 * p.ApexY)
//                    |> Array.unzip
//                match caruanaAlgorithm xData yData with 
//                | Some caruanaEst -> 
//                    let skew = skewBy caruanaEst.MeanX  (caruanaEst.STD**2.) p.XData p.YData
//                    let tau  = estTau caruanaEst.STD skew
//                    Some (createEstimatedMoments p.ApexY caruanaEst.MeanX (caruanaEst.STD**2.) caruanaEst.STD skew tau)
//                | Option.None            -> 
//                    methodOfMoments p 
      
 
//    module HULQ = 

//        open Integration
//        open ParameterEstimation
//        open LevenbergMarquardt

//        type PeakModel = 
//            | Gaussian of Model
//            | EMG of Model 

//        type FittedPeak = {
//            Model                       : PeakModel 
//            EstimatedParams             : float []
//            StandardErrorOfPrediction   : float
//            YPredicted                  : float []
//            }

//        ///
//        let createFittedPeak model estimatedParams standardErrorOfPrediction yPredicted = {
//            Model                       = model                       
//            EstimatedParams             = estimatedParams            
//            StandardErrorOfPrediction   = standardErrorOfPrediction  
//            YPredicted                  = yPredicted                 
//            }

//        type QuantifiedPeak = {
//            Model                       : PeakModel option 
//            YPredicted                  : float []
//            EstimatedParams             : float []
//            StandardErrorOfPrediction   : float
//            Area                        : float
//            MeasuredApexIntensity       : float
//            }

//        ///
//        let createQuantifiedPeak model yPredicted estimatedParams standardErrorOfPrediction area measuredApexIntensity = {
//            Model                       = model                       
//            YPredicted                  = yPredicted                  
//            EstimatedParams             = estimatedParams             
//            StandardErrorOfPrediction   = standardErrorOfPrediction   
//            Area                        = area      
//            MeasuredApexIntensity       = measuredApexIntensity
//            }

//        /// Return Option
//        let getPeakBy x (peaks:SignalEntity<'a> []) =
//            let isPartOfPeak x (p:SignalEntity<'a>)  = p.Start < x && p.End > x        
//            match Array.tryFind (isPartOfPeak x) peaks with
//            | Some p        ->  p
//            | Option.None   -> 
//                Array.minBy (fun p -> abs(0.5*(p.Start + p.End) - x)) peaks 
            
//        ///
//        let tryFit model solverOptions xData yData =
//            let modelF = 
//                match model with 
//                | Gaussian f -> f
//                | EMG f      -> f 
                
//            try
//                let estParams = estimatedParamsVerbose modelF solverOptions 0.001 10.0 xData yData |> Seq.last                                  
//                let y' = Array.map (fun xValue -> modelF.GetFunctionValue estParams xValue) xData
//                let sEoE_Gauss = standardErrorOfPrediction (float solverOptions.InitialParamGuess.Length) y' yData
//                Some (createFittedPeak model (estParams.ToArray()) sEoE_Gauss y')
//            with 
//            | ex -> 
//                printfn "%A" ex
//                Option.None      

//        ///
//        let tryFitGaussian initAmp initMeanX initStdev xData yData =
//            let solverOptions = createSolverOption 0.00001 0.001 1000 [|initAmp;initMeanX;initStdev|]
//            tryFit (Gaussian Table.gaussModel) solverOptions xData yData

//        ///
//        let tryFitEMG initAmp initMeanX initStdev initTau xData yData =
//            let solverOptions = createSolverOption 0.00001 0.001 1000 [|initAmp;initMeanX;initStdev;initTau|]
//            tryFit (EMG Table.emgModel) solverOptions xData yData

//        ///
//        let selectModel (models: FittedPeak []) = 
//            Array.minBy (fun (x:FittedPeak) -> x.StandardErrorOfPrediction) models

//        ///
//        let integralOfGaussian (coeffs:float []) =
//            sqrt(2.) * coeffs.[0] * coeffs.[2] * sqrt(Math.PI)

//        ///
//        let integralOfEMGBy meanX sigma tau intensityAtx xPos =    
//            let fst = intensityAtx * (sqrt(Math.PI)*tau)
//            let intg =         
//                let firstT = sqrt(Math.PI) / 2.
//                let errorF = 
//                    ((xPos-meanX)/(sqrt(2.)*sigma)) - (sigma/(sqrt(2.)*tau))
//                    |> FSharp.Stats.SpecialFunctions.Errorfunction.Erf
//                firstT * (1.+errorF)
//            let exp = exp ( 0.5*((sigma/tau)**2.) - ((xPos-meanX)/tau) )
//            fst / (intg * exp)

//        ///
//        let calcArea (fit: FittedPeak) =
//            ///
//            let integralOfGaussianOf (fit:FittedPeak) =
//                integralOfGaussian fit.EstimatedParams
//            ///
//            let integralOfEMGOf (fit:FittedPeak) =
//                let pF = 
//                    match fit.Model with 
//                    | EMG f -> f
//                    | _     -> failwith "fittedPeak was not fitted with emg"
//                let intensityAtx = pF.GetFunctionValue (fit.EstimatedParams |> vector) fit.EstimatedParams.[1]
//                integralOfEMGBy fit.EstimatedParams.[1] fit.EstimatedParams.[2] fit.EstimatedParams.[3] intensityAtx fit.EstimatedParams.[1]
//            match fit.Model with 
//            | Gaussian _ -> integralOfGaussianOf fit
//            | EMG      _ -> integralOfEMGOf fit

//        module SecondDerivative = 

//            ///
//            let toSignalEntity (p:IdentifiedPeak ) :SignalEntity<'a option> = 
//                let peakToQuantify = {
//                    ApexY = p.Apex.YVal 
//                    ApexX = p.Apex.XVal
//                    XData = p.XData
//                    YData = p.YData
//                    ParameterEstimates = Option.None
//                    }

//                {
//                    Start                   = p.LeftEnd.XVal 
//                    End                     = p.RightEnd.XVal
//                    XData                   = p.XData
//                    YData                   = p.YData
//                    PeaksToQuantify         = [peakToQuantify]
//                    HasParameterEstimates   = false
//                    HasMultiplePeaks        = false
//                }

//        module Wavelet = 

//            ///
//            let toSignalEntity (p:PeakGroup) :SignalEntity<(EstimatedMoments)> = 
//                let peakToQuantify (g:Gaussian) : PeakToQuantify<EstimatedMoments> = 
//                    let x,y = g.Trace |> Array.unzip
//                    let tau = 
//                    let moments = createEstimatedMoments g.Amplitude g.XLoc g.Stdev (sqrt g.Stdev)
//                    {
//                    ApexY = g.XLoc 
//                    ApexX = g.Amplitude
//                    XData = x
//                    YData = y
//                    ParameterEstimates = (Option.Some moments)
//                    }

//                let x,y = p.Data |> Array.unzip 
//                {
//                    Start                   = p.Start 
//                    End                     = p.End
//                    XData                   = x
//                    YData                   = y
//                    PeaksToQuantify         = (p.Fits |> List.map peakToQuantify) 
//                    HasParameterEstimates   = true
//                    HasMultiplePeaks        = p.Convolved
//                }

//        ///
//        let quantifySinglePeak (getMoments:PeakToQuantify<'a> -> EstimatedMoments option) (p: PeakToQuantify<'a>) =
//            match getMoments p with 
//            | Some moments -> 
//                let gaussian = tryFitGaussian moments.ModeY moments.MeanX moments.Std  p.XData p.YData
//                let emg      = tryFitEMG moments.ModeY moments.MeanX moments.Std moments.Tau p.XData p.YData
//                match gaussian, emg with 
//                | Some g, Some emg -> 
//                    let finalFit   = selectModel [|g;emg|]
//                    let area       = calcArea finalFit
//                    createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction area p.ApexY
//                | Some finalFit, Option.None  | Option.None, Some finalFit -> 
//                    let area       = calcArea finalFit
//                    createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction area p.ApexY
//                | _                       ->
//                    let area = trapezEstAreaOf p.XData p.YData
//                    createQuantifiedPeak Option.None [||] [||] nan area p.ApexY
//            | Option.None   -> 
//                let area = trapezEstAreaOf p.XData p.YData
//                createQuantifiedPeak Option.None [||] [||] nan area p.ApexY         

//        ///
//        let quantifyConvolvedPeak (getMoments:SignalEntity<'a> -> (EstimatedMoments option) list) (p: SignalEntity<'a>) =
//            let moments = getMoments p |> List.choose id
//            match moments with 
//            | [] ->
//                let area = trapezEstAreaOf p.XData p.YData
//                createQuantifiedPeak Option.None [||] [||] nan area (p.YData |> Array.max)         
//            | _ -> 
//                let area = trapezEstAreaOf p.XData p.YData
//                createQuantifiedPeak Option.None [||] [||] nan area (p.YData |> Array.max)       
//            //match getMoments p with 
//            //| Some moments -> 
//            //    let gaussian = tryFitGaussian moments.ModeY moments.MeanX moments.Std  p.XData p.YData
//            //    let emg      = tryFitEMG moments.ModeY moments.MeanX moments.Std moments.Tau p.XData p.YData
//            //    match gaussian, emg with 
//            //    | Some g, Some emg -> 
//            //        let finalFit   = selectModel [|g;emg|]
//            //        let area       = calcArea finalFit
//            //        createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction area p.ApexY
//            //    | Some finalFit, Option.None  | Option.None, Some finalFit -> 
//            //        let area       = calcArea finalFit
//            //        createQuantifiedPeak  (Some finalFit.Model) finalFit.YPredicted finalFit.EstimatedParams finalFit.StandardErrorOfPrediction area p.ApexY
//            //    | _                       ->
//            //        let area = trapezEstAreaOf p.XData p.YData
//            //        createQuantifiedPeak Option.None [||] [||] nan area p.ApexY
//            //| Option.None   -> 
//            //    let area = trapezEstAreaOf p.XData p.YData
//            //    createQuantifiedPeak Option.None [||] [||] nan area p.ApexY         
                
//        /////
//        //let quantifyBy (getMoments:PeakToQuantify<'a> -> EstimatedMoments option) (getMoments:SignalEntity<'a> -> EstimatedMoments option) x (peaks: SignalEntity<'a> [] ) =
//        //    let peakToQuantify = getPeakBy x peaks  
//        //    match peakToQuantify.HasMultiplePeaks with
//        //    | false  -> 
//        //        let p = peakToQuantify.PeaksToQuantify |> Seq.head
//        //        quantifyPeak getMoments p
//        //    | true-> 
//        //        let p = peakToQuantify.PeaksToQuantify |> Seq.head
//        //        quantifyPeak p

            

<|start_filename|>src/BioFSharp.Mz/Peak.fs<|end_filename|>
﻿namespace BioFSharp.Mz

type IPeak =    
    abstract member Mz        : float
    abstract member Intensity : float   


[<Struct>]
type Peak(mz:float,intensity:float) = 
        member this.Mz = mz
        member this.Intensity = intensity
        
        interface IPeak with
            member this.Mz = mz
            member this.Intensity = intensity

type PeakFamily<'a> = {
    MainPeak       : 'a
    DependentPeaks : 'a list   
}

module Ions = 
    open System.Linq

    [<System.Flags>]
    type IonTypeFlag =
        | Precursor     =  2
        | A             =  4
        | B             =  8
        | C             =  16
        | X             =  32
        | Y             =  64
        | Z             =  128
        | lossH2O       =  256
        | lossNH3       =  512
        | Immonium      = 1024
        | Neutral       = 2048
        | Diagnostic    = 4096
        | Unknown       = 8192
      
    let createIonTypeList (ionTypeFlag: IonTypeFlag) =
        System.Enum.GetValues(ionTypeFlag.GetType()).Cast<IonTypeFlag>().Where(ionTypeFlag.HasFlag)

    let hasFlag (ionType:IonTypeFlag) (ionType':IonTypeFlag) = 
        (int ionType &&& int ionType') <> 0


module Peaks =
 
    ///
    let createPeak mzData intensityData = 
        let pk = Peak(mzData, intensityData)
        pk
    

    ///
    let createPeakFamily mainPeak dependentPeaks = {
        MainPeak       = mainPeak
        DependentPeaks = dependentPeaks  
        }

module TaggedMass = 
    
    [<Struct>]
    type TaggedMass (iontype:Ions.IonTypeFlag,mass:float) = 
            member this.Iontype = iontype
            member this.Mass = mass
            
    ///
    let createTaggedMass iontype mass = 
        let tM = TaggedMass(iontype,mass)
        tM   
    
    ///
    let createTaggedH2OLoss iontype mass =
        let tM = TaggedMass(iontype+Ions.IonTypeFlag.lossH2O,mass)
        tM   

    ///
    let createTaggedNH3Loss iontype mass =
        let tM = TaggedMass(iontype+Ions.IonTypeFlag.lossNH3,mass)
        tM   

    
module TaggedPeak =
    
    [<Struct>]
    type TaggedPeak (iontype:Ions.IonTypeFlag,mz:float,intensity:float) = 
            member this.Iontype = iontype
            member this.Mz = mz
            member this.Intensity = intensity
            interface IPeak with
                member this.Mz = mz
                member this.Intensity = intensity
            
    ///
    let createTaggedPeak iontype mzData intensityData = 
        let tPk = TaggedPeak(iontype,mzData,intensityData)
        tPk    
    
    /// 
    let createTaggedPeakOf (taggedMass: TaggedMass.TaggedMass) (intensityPredF: Ions.IonTypeFlag -> float) =
        createTaggedPeak taggedMass.Iontype taggedMass.Mass (intensityPredF taggedMass.Iontype)


<|start_filename|>src/BioFSharp.Mz/SignalDetection.fs<|end_filename|>
﻿namespace BioFSharp.Mz
// The module contains a modified version of a wavelet-based peak picking approach
// published by French et. al, available via doi: 10.1021/pr500886y.
// The original work was published using the Apache v2 license https://www.apache.org/licenses/LICENSE-2.0

open BioFSharp
open BioFSharp.IO

open System
open FSharpAux
open AminoAcids 
open ModificationInfo

open FSharp.Stats

module SignalDetection =
    
    open BioFSharp.Mz.PeakArray

    module Care =    

        ///
        let getColLowBound  (xData: float[]) centerIdx mzTol = 
            let rec loop  (xData: float []) center i mzTol =
                if i <= 0 then 0
                else                                              
                    match (xData.[center] - xData.[i]) with 
                    | x when x >= mzTol -> i+1   
                    | x when x <= mzTol -> loop xData center (i-1) mzTol
            loop xData centerIdx (centerIdx-1) mzTol

        /// 
        let getColHighBound (xData: float []) centerIdx mzTol = 
            let rec loop (xData: float []) center i mzTol = 
                if i >= xData.Length-1 then xData.Length-1
                else                                             
                    match xData.[i] - xData.[center] with 
                    | x when x >= mzTol -> i-1
                    | x when x <= mzTol -> loop xData center (i+1) mzTol
            loop xData centerIdx (centerIdx+1) mzTol
    
        ///
        let calc targetX state currentX currentY =
            let m,i = state
            (m+currentX*currentY,i+currentY)


        /// Returns a Array containing the distances between adjacent mz values of the input array.
        let createXspacing (mzData:float[]) =
            let xSpacing = Array.create (mzData.Length) 0.
            for mzi=1 to mzData.Length-1 do
                    xSpacing.[mzi] <- mzData.[mzi] - mzData.[mzi-1]
            xSpacing.[0] <- xSpacing.[1]
            xSpacing.[mzData.Length-1] <- xSpacing.[mzData.Length-2]
            xSpacing

        /// Helper function to accumulate Array values
        let accumulate fromIdx toIdx startV f (data:'T[]) =
            let rec loop i acc =
                if i = toIdx then 
                    acc
                else
                    loop (i+1) (f acc data.[i])
            loop fromIdx startV  

        /// Helper function for determining the score in a (sorted) vector at a given percentile
        /// end of getScoreAtPercentile. Allow passing of length of vector in case you only 
        /// want a slice of the first portion of a vector;
        /// perc should not be a fraction (e.g. 5th per centile = 5.0)
        let scoreAtPercentile (perc:float) nTot (sortedData: float [])  = 
            if nTot = 0 then 0.
            else
            let nBelow = (float (nTot-1)) * perc / 100. //Anzahl mal percentile 
            if (ceil nBelow) = nBelow 
                then sortedData.[int nBelow] //whole number
            else 
                let lowIndex = floor nBelow
                let highIndex = ceil nBelow
                let fraction = nBelow - lowIndex
                let low = sortedData.[int lowIndex]
                let high = sortedData.[int highIndex]
                (low + (high - low) * fraction)
    
    // TODO: add to FSharp.Stats.Signal
    module Padding = 
        
        ///
        type PaddingParameters = { 
            PaddingYValue           : float
            MaximumPaddingPoints    : int option
            MzTolerance             : float 
            WindowSize              : int
            SpacingPerc             : float 
            }
           
        ///   
        let createPaddingParameters paddingYValue maximumPaddingPoints mzTolerance windowSize spacingPerc = {
            PaddingYValue           = paddingYValue           
            MaximumPaddingPoints    = maximumPaddingPoints    
            MzTolerance             = mzTolerance             
            WindowSize              = windowSize              
            SpacingPerc             = spacingPerc             
            }      

        ///
        let paddDataWith paddYValue maxPaddingPoints (mzTol: float) (windowSize: int) (spacing_perc:float) (xData: float []) (yData: float []) = 
            //////////////////////////////////////////////////////////////////////////
            //Step 1. Estimate the spacing between adjacent mzPoints (mz- (mz-1))
            let xSpacing = Care.createXspacing xData
            //////////////////////////////////////////////////////////////////////////
            //Step 2. Estimate the spacing threshold in bins
            let mutable windowSize = windowSize
            if (windowSize > xSpacing.Length-1) then
                windowSize <- (xSpacing.Length-1) / 2
            let hf_window = windowSize / 2
            windowSize <- 2 * hf_window
            let nSpacingBins = ((xSpacing.Length-1) / (windowSize + 1))
            let spacings = Array.zeroCreate (nSpacingBins+1) 
            for i = 0 to nSpacingBins-1 do
                let windowLow = i * windowSize
                let mutable windowHigh = windowLow + windowSize //not inclusive
                if i = nSpacingBins-1 
                    then windowHigh <- xSpacing.Length
                let nTot = (windowHigh - windowLow)
                let unsortedData = Array.zeroCreate nTot //nTot ist = windowsize, es sei den die Distanz von windowlow zum array Ende ist größer als windowsize
                for j = windowLow to windowHigh-1 do
                    unsortedData.[j-windowLow] <- xSpacing.[j] 
                let sortedData = unsortedData |> Array.sortDescending 
                let mutable spacingAtPerc = Care.scoreAtPercentile spacing_perc nTot sortedData  //calcs Noisethreshold in bin i                 
                spacings.[i] <- spacingAtPerc
            //////////////////////////////////////////////////////////////////////////
            //Step 3. If the mzSpacing to the average spacing in the selected bin is greater than 1. add mz points with zero intensity to the data
            let maxSpacingToAvgSpacing = 1.
            let xyData = Array.zip xData yData 
            let paddedData = ResizeArray<float*float>()
            for mzi = 1 to xyData.Length-2 do
                let spacingBin =  //TODO zu care
                    if  nSpacingBins = 0 then 0
                    else
                        let tmp = mzi/windowSize
                        match tmp with
                        | tmp when tmp > nSpacingBins-1 -> nSpacingBins-1 
                        | _ -> tmp
                let spacingToAvgSpacing = xSpacing.[mzi+1] / spacings.[spacingBin]                
                if spacingToAvgSpacing > maxSpacingToAvgSpacing && xSpacing.[mzi+1] > mzTol then                    
                        paddedData.Add xyData.[mzi]
                        let maxPossible = (int (xSpacing.[mzi+1] / mzTol) - 1 )
                        let maxPaddings = 
                            match maxPaddingPoints with 
                            | None     -> maxPossible
                            | Some max -> 
                                if max >= maxPossible then maxPossible
                                else max                      
                        for i=1 to maxPaddings do
                            paddedData.Add (xData.[mzi] + ((float i)*mzTol), paddYValue)
                        let leftMostMz = fst paddedData.[paddedData.Count-1]
                        for i=1 to maxPaddings do
                            let newValue = fst (xData.[mzi+1] - ((float i)*mzTol), paddYValue)
                            if newValue > leftMostMz then 
                                paddedData.Add (xData.[mzi+1] - ((float i)*mzTol), paddYValue)                    
                else
                    paddedData.Add xyData.[mzi]
            let nTotFinal = paddedData.Count
            let xDataPadded = Array.zeroCreate nTotFinal
            let yDataPadded = Array.zeroCreate nTotFinal            
            paddedData 
            |> Seq.iteri (fun i (x,y) -> 
                            xDataPadded.[i] <- x
                            yDataPadded.[i] <- y
                        )  
            xDataPadded, yDataPadded   

        ///
        let paddDataExtensively paddValue (mzTol: float) (windowSize: int) (spacing_perc:float) (xData: float []) (yData: float []) =             
            paddDataWith paddValue None mzTol windowSize spacing_perc xData yData
        
        ///
        let paddDataBy (paddingParameters:PaddingParameters) (xData: float []) (yData: float []) = 
            paddDataWith paddingParameters.PaddingYValue paddingParameters.MaximumPaddingPoints paddingParameters.MzTolerance 
                paddingParameters.WindowSize paddingParameters.SpacingPerc xData yData

    module Wavelet =

        [<Struct>]
        type RidgeLine = {
            Row: int
            Col: int
            }
        
        ///
        let createRidgeLine row col = {
            Row=row;
            Col=col; 
            }

        ///
        type WaveletParameters = {
            NumberOfScales        : int
            YThreshold            : float
            MzTolerance           : float
            SNRS_Percentile       : float
            MinSNR                : float 
            RefineMZ              : bool
            SumIntensities        : bool
            }

        //toCentroid ricker2d nScales yYhreshold mzTol (*spacing_perc*) snrs_perc minSnr mzData intensityData
        ///
        let createWaveletParameters numberOfScales yThreshold mzTolerance sNRS_Percentile minSNR refineMZ sumIntensities = {
            NumberOfScales  = numberOfScales
            YThreshold      = yThreshold 
            MzTolerance     = mzTolerance
            SNRS_Percentile = sNRS_Percentile
            MinSNR          = minSNR
            RefineMZ        = refineMZ
            SumIntensities  = sumIntensities
            }

        /// Helperfunction to calculate the mz values left and right from the target value which will be included in the computation of 
        /// the CWT correlation coefficient
        let getNpointsX centerIdx scaling maxMZwindow (nPoints: int [,,]) (mzData:float[])  =
            ///
            let rec getNPointsLeft acc counter =  
                if counter >= 0 then
                    if (mzData.[centerIdx] - mzData.[counter]) > maxMZwindow  then 
                         nPoints.[0, scaling, centerIdx] <- acc                     
                    else getNPointsLeft (acc+1) (counter-1) 
                else nPoints.[0, scaling, centerIdx] <- acc          
            getNPointsLeft 0 centerIdx 
            ///
            let rec getNPointsRight acc' counter  =
                if counter < mzData.Length then      
                    if ((mzData.[counter] - mzData.[centerIdx]) > maxMZwindow) then       
                        nPoints.[1, scaling, centerIdx] <- acc'
                    else getNPointsRight (acc'+1) (counter+1) 
                else 
                     nPoints.[1, scaling, centerIdx] <- acc'
            getNPointsRight 0 centerIdx   
                
        ///
        let createPadding paddingPoints (mzDataPadded:float[]) (intensityDataPadded:float[]) (mzData:float[]) (intensityData: float[]) = 
            for i = paddingPoints to mzData.Length - 1 + paddingPoints do 
                mzDataPadded.[i] <- mzData.[i-paddingPoints]
                intensityDataPadded.[i] <- intensityData.[i-paddingPoints]

        /// Helperfunction to get mz back of Columnnumber
        let convertColToMz (mzData: float []) (col: int) =
            let mapIndex = col / 2
            if ( col % 2 = 1)
                then (mzData.[mapIndex] + mzData.[mapIndex+1]) / 2.0
            else 
                mzData.[mapIndex] 

        ///
        let getScales nScales (scalings:float []) (*spacing_perc*) windowMaxLow windowMaxHigh (mzData:float[])  (intensityData: float[]) (xSpacingAvg:float[]) (nPoints: int [,,]) (xspaced:float[]) =
            for mzi = 1 to mzData.Length-2 do
                let windowLow' = if mzi < windowMaxLow then 0 else mzi - windowMaxLow
                let windowHigh' = if mzi >= mzData.Length - windowMaxHigh then mzData.Length - 1 else mzi + windowMaxHigh
                let nTot = windowHigh' - windowLow' |> float
                xSpacingAvg.[mzi] <- ((Care.accumulate windowLow' windowHigh' 0. (+) xspaced) / nTot)
            for mzi = 1 to mzData.Length-2 do
                let mutable scalesToInclude = nScales    
                if intensityData.[mzi] = 1. ||intensityData.[mzi] < (0.75*intensityData.[mzi-1]) || intensityData.[mzi] < 0.75*intensityData.[mzi+1] then 
                        scalesToInclude <- 1
                // figure out the number of wavelet points you'll need to sample for each m/z point
                for i = 0 to scalesToInclude-1 do 
                    let maxMZwindow = 
                        let tmp = xSpacingAvg.[mzi] * scalings.[i] *  3.0 
                        if tmp > 3. then 3. else tmp 
                    getNpointsX mzi i maxMZwindow nPoints mzData

        /////
        //let ricker2d (mzDataPadded:float []) (waveletData:float []) (centerIdx:int) (nPointsLeft:int) (nPointsRight:int) (focusedPaddedMzValue:float) (width:float) =  //(PadMz, paddedCol, nPointsLeft, nPointsRight, param1, param2, PadMz[paddedCol], waveletData)
        //    let computation cnt i =
        //        let vec = mzDataPadded.[i]-focusedPaddedMzValue
        //        let tsq = vec*vec
        //        let modi = 1. - tsq / (width * width)
        //        let gauss = exp(-1. * tsq / (2.0 *(width * width)) )
        //        waveletData.[cnt] <- (  2.0 / ( sqrt(3.0 * width) * (sqrt(sqrt(3.141519))) ) ) * modi * gauss
        //    let rec ricker cnt i = 
        //        if i = centerIdx+nPointsRight then 
        //             computation cnt i
        //        else computation cnt i
        //             ricker (cnt+1) (i+1)
        //    if nPointsLeft - nPointsRight >= waveletData.Length 
        //        then printfn "invalid input Parameters for ricker2d"
        //    ricker 0 (centerIdx-nPointsLeft)
        //    waveletData

        ///
        let ricker2d (mzDataPadded:float []) (waveletData:float []) (centerIdx:int) (nPointsLeft:int) (nPointsRight:int) (focusedPaddedMzValue:float) (width:float) =  //(PadMz, paddedCol, nPointsLeft, nPointsRight, param1, param2, PadMz[paddedCol], waveletData)
            let inline computation width2 param2 param3 cnt i =
                let vec = mzDataPadded.[i]-focusedPaddedMzValue
                let tsq = vec*vec
                let modi = 1. - tsq / width2
                let gauss = exp(-1. * tsq / param2 )
                waveletData.[cnt] <- param3 * modi * gauss
            let rec ricker width2 param2 param3 cnt i = 
                if i = centerIdx+nPointsRight then 
                        computation width2 param2 param3 cnt i
                else computation width2 param2 param3 cnt i
                     ricker width2 param2 param3 (cnt+1) (i+1)
            if nPointsLeft - nPointsRight >= waveletData.Length 
                then printfn "invalid input Parameters for ricker2d"
            let width2 = (width * width)
            let param2 = (2.0 * width2)
            let param3 = 2.0 / ( sqrt(3.0 * width) * (sqrt(sqrt(3.141519))))
            ricker width2 param2 param3 0 (centerIdx-nPointsLeft)
            waveletData 
        
        ///       
        let inline calcCorrWaveletAndDataMatrix (waveletF: float [] -> float [] -> int -> int -> int -> float -> float -> float [])  nScales paddingPoints (mzDataPadded:float[]) (intensityDataPadded:float[]) (mzData:float[]) (intensityData: float[])(* (waveletData:float[])*) (nPoints: int[,,]) (xSpacingAvg: float []) (scalings: float []) (corrMatrix: float [,]) = 
            let calcCorrAtScale i =               
                let waveletData' = Array.create (paddingPoints) 0. 
                let currentscaling = scalings.[i]
                let mutable matrixIndex = 2 
                let calcCorrelationAtPosX targetCol =
                    let nPointsLeft = nPoints.[0,i,targetCol]
                    let nPointsRight = nPoints.[1,i,targetCol]
                    if nPointsLeft = 0 && nPointsRight = 0 then 
                        matrixIndex <- (matrixIndex+2)
                    else
                        let centerPaddedColIdx = targetCol + paddingPoints 
                        let width = xSpacingAvg.[targetCol] * currentscaling
                        // ricker, schreibt in "waveletdata", ein Array dem die Abstände zwischen den punkten zugrunde liegen. 
                        // und gibt die Werte der Waveletfunktion in Abhängigkeit des fokussierten m/z values wieder 
                        let correlation = waveletF mzDataPadded waveletData' centerPaddedColIdx nPointsLeft nPointsRight mzDataPadded.[centerPaddedColIdx]  width
                        //daraufhin wird die Korrelationsmatrix gefüllt unzwar an der Spalte des aktuellen scalings gepaart mit einem aufsteigenden Index welcher der Nr des m/z values
                        //im Urpsprungs m/z Array entspricht. Dafür wird der "waveletData" Array durchlaufen und mit der Intensität an der entsprechenden Stelle multipliziert, die Summe der Ergebnisse wird in 
                        //die Korrelationsmatrix eingetragen. Es handelt sich um den CWT Koeffizienten
                        let startpoint = centerPaddedColIdx - nPointsLeft                          
                        let mutable acc = 0.
                        for k = 0 to nPointsLeft+nPointsRight do
                            acc <- acc + (correlation.[k] * intensityDataPadded.[k + startpoint] ) 
                            //corrMatrix.[i, matrixIndex] <- (corrMatrix.[i, matrixIndex] + (correlation.[k] * intensityDataPadded.[k + startpoint] ) )          
                        corrMatrix.[i, matrixIndex] <- acc
                        matrixIndex <- matrixIndex+1

                        if targetCol < mzData.Length-1 then //falls j = mzData.Length.1 dann ist man am ende des m/z Arrays angekommen und die corrMatrix ist fertig            // falls nicht : Dieser Vorgang wird wiederholt unzwar wird dazu ein waveletdataarray erzeugt, der auf den Punkten zwischen den m/z werden basiert.
                            let moverzShift =  ( (mzDataPadded.[centerPaddedColIdx] + mzDataPadded.[centerPaddedColIdx+1]) / 2.0 )
                            let correlation = waveletF mzDataPadded waveletData' centerPaddedColIdx nPointsLeft nPointsRight moverzShift width                                
                            let mutable acc = 0.
                            for k = 0 to nPointsLeft+nPointsRight do                
                                acc <- acc + (correlation.[k] * intensityDataPadded.[k + startpoint] ) 
                                //corrMatrix.[i, matrixIndex] <- (corrMatrix.[i, matrixIndex] + (correlation.[k] * intensityDataPadded.[k + startpoint] ) )
                            corrMatrix.[i, matrixIndex] <- acc
                            matrixIndex <- matrixIndex+1 
                let calcCorrAt j =
                    if i>0 then 
                        if intensityData.[j] < 0.75*intensityData.[j-1] || intensityData.[j] < 0.75*intensityData.[j+1] then 
                                 matrixIndex <- (matrixIndex+2) 
                            else calcCorrelationAtPosX j                   
                    else
                        calcCorrelationAtPosX j 
                for j = 1 to mzData.Length-2 do
                   calcCorrAt j

            [0.. nScales-1]
            |> PSeq.iter calcCorrAtScale
            //[for i in [0.. nScales-1] -> 
            //    async {return calcCorrAtScale i }]
            //|> Async.Parallel
            //|> Async.RunSynchronously

        ///
        let getSNRFilteredPeakLines noise_perc minSnr mzTol nScales mzData (allLines:Collections.Generic.List<RidgeLine>) (corrMatrix: float[,]) =                                           
            let corrMatrixLength = corrMatrix.[0,*].Length
            /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
            // step 1: find maxima in each column (row = scales, columns = m/z)
            /// contains scaling with the highest correlationvalue for each Column
            let maxRowofCol = Array.create (corrMatrixLength) 0 
            for i = 0 to corrMatrixLength-1 do             
                let mutable corrMax = 0.0            
                for j = 0 to nScales-1 do 
                    if corrMatrix.[j,i] > corrMax then
                        corrMax <- corrMatrix.[j,i]
                        maxRowofCol.[i] <- j //Eintrag des scalings mit der höchsten correlation
            /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
            // step 2, setup bins of 300 points and calculate noise threshold within each bin
            let mutable windowSize = 300
            if (windowSize > corrMatrixLength) 
                then windowSize <- corrMatrixLength / 2
            let hf_window = windowSize / 2
            windowSize <- 2 * hf_window
            let nNoiseBins = corrMatrixLength / windowSize + 1 //number of NoiseBins
            let noises = Array.create nNoiseBins 0.0
            for i = 0 to nNoiseBins-1 do
                let windowLow = i * windowSize  // inclusive
                let mutable windowHigh = windowLow + windowSize //not inclusive
                
                if i = nNoiseBins-1 
                    then windowHigh <- corrMatrixLength
                let nTot = windowHigh - windowLow                 
                let unsortedData = Array.create nTot 0.0 //nTot ist = windowsize, es sei den die Distanz von windowlow zum array Ende ist größer als windowsize
                for j = windowLow to windowHigh - 1 do
                    unsortedData.[j-windowLow] <- corrMatrix.[0, j] // first row of correlation matrix                
                let sortedData = unsortedData |> Array.sort
                let mutable noise = Care.scoreAtPercentile noise_perc nTot sortedData  //calcs Noisethreshold in bin i 
                if noise < 1.0 then noise <- 1.0
                noises.[i] <- noise
            let interpolatedXpoints = Array.create corrMatrixLength 0.0 
            for i = 0 to corrMatrixLength-1 do
                let interpolatedpointOfI = convertColToMz mzData i //mz Value of Colum
                interpolatedXpoints.[i] <- interpolatedpointOfI
            /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
            // step 3, find local maxima that are separated by at least mzTol_
            let findLocalMaxima i = 
                // get the indices for the lower and upper bounds   
                let lowBound  = Care.getColLowBound interpolatedXpoints i (mzTol/2.)
                let highBound = Care.getColHighBound interpolatedXpoints i (mzTol/2.)
                // Now all m/z Values in the space between lowBound and highBound get analyzed and 
                // the Column of the m/z with the  maxCorr is identified
                let mutable maxCorr = 0.0
                let mutable maxCol = 0
                for j = lowBound to highBound do 
                    let row = maxRowofCol.[j]
                    if corrMatrix.[row,j] > maxCorr then
                        maxCorr <- corrMatrix.[row,j]
                        maxCol <- j              
                let noiseBin = 
                    let tmp = maxCol/windowSize
                    match tmp with
                    | tmp when tmp > nNoiseBins-1 -> nNoiseBins-1 
                    | _ -> tmp    
                let snr = maxCorr / noises.[noiseBin]
                let row = maxRowofCol.[maxCol]
                if snr > minSnr && row > 0 then // 
                    let nLines = allLines.Count
                    if ( nLines > 0) then
                        let mzNewLine = convertColToMz mzData maxCol 
                        let mzPrevLine = convertColToMz mzData allLines.[(nLines-1)].Col
                        let mzDiff = mzNewLine - mzPrevLine
                        let corrPrev = corrMatrix.[allLines.[nLines-1].Row, allLines.[nLines-1].Col]
     
                        if mzDiff > mzTol then 
                            let newLine = createRidgeLine maxRowofCol.[maxCol] maxCol
                            allLines.Add newLine
                        elif maxCorr > corrPrev then 
                            let newLine = createRidgeLine maxRowofCol.[maxCol] maxCol
                            allLines.[allLines.Count-1] <- newLine
                    else 
                        let newLine = createRidgeLine maxRowofCol.[maxCol] maxCol
                        allLines.Add newLine
            let rec findLocalMaximaLoop i (corrMatrix: float [,]) = 
                if i = corrMatrixLength-3 then 
                    findLocalMaxima i
                else                  
                        findLocalMaxima i
                        let correlationVal = corrMatrix.[maxRowofCol.[i],i]
                        if  correlationVal < corrMatrix.[maxRowofCol.[i-1],i-1] ||
                            correlationVal < corrMatrix.[maxRowofCol.[i-2],i-2] ||
                            correlationVal < corrMatrix.[maxRowofCol.[i+1],i+1] ||
                            correlationVal < corrMatrix.[maxRowofCol.[i+2],i+2] then 
                                findLocalMaximaLoop (i+1) corrMatrix
                        else 
                            findLocalMaxima i
                            findLocalMaximaLoop (i+1) corrMatrix
            findLocalMaximaLoop 2 corrMatrix 

        /// 
        let refinePeaks refineMz sumIntensities yThreshold mzTol (scalings: float[]) (mzData:float []) (intensityData:float [])  (allLines: ResizeArray<RidgeLine>) (xSpacingAvg:float []) =          
            let xData = ResizeArray<float>(2500)
            let yData = ResizeArray<float>(2500)
            if allLines.Count = 0 then 
                let finalX = xData.ToArray()
                let finalY = yData.ToArray()
                finalX, finalY 
            else
                for i = 0 to allLines.Count-1 do
                    let mzColIdx = allLines.[i].Col / 2
                    let row = allLines.[i].Row
                    let currentScaling = scalings.[row]
                    let offset = 
                        let offsetTmp = currentScaling * xSpacingAvg.[int (allLines.[i].Col / 2)]
                        if offsetTmp > mzTol then
                            mzTol
                        else 
                            offsetTmp 
                    let startFittingPoint = Care.getColLowBound mzData mzColIdx offset                   
                    let endFittingPoint = Care.getColHighBound mzData mzColIdx offset
                    // sum up the intensity and find the highest point. If there are multiple
                    // points with the maxIntensity value, take the one with the highest m/z.
                    let mutable maxIntensityMZ = 0.0
                    let mutable maxIntensity = 0.0 
                    let mutable weightedMzSum = 0.0
                    let mutable intensitySum = 0.0 
                    for j = startFittingPoint to endFittingPoint do 
                        intensitySum <- intensitySum + intensityData.[j] 
                        weightedMzSum <- weightedMzSum + (mzData.[j] * intensityData.[j])
                        if intensityData.[j] >= maxIntensity
                            then
                                maxIntensity   <- intensityData.[j]
                                maxIntensityMZ <- mzData.[j] 
                    let refCentLength = xData.Count    
                    let currentMz = 
                        if refineMz then 
                            let weightedMeanMz = weightedMzSum / intensitySum
                            weightedMeanMz 
                        else
                            maxIntensityMZ
                    let currentIntensity =
                        if sumIntensities then 
                            intensitySum
                        else
                            maxIntensity
                    if refCentLength = 0 && maxIntensity > yThreshold then
                        xData.Add currentMz
                        yData.Add currentIntensity
                    elif refCentLength > 0 && (currentMz - xData.[refCentLength-1]) < mzTol && (currentIntensity > yData.[refCentLength-1] && maxIntensity > yThreshold ) then 
                        xData.[refCentLength-1] <- currentMz
                        yData.[refCentLength-1] <- currentIntensity             
                    elif refCentLength > 0 && currentMz <> xData.[refCentLength-1] && (currentMz - xData.[refCentLength-1] ) > mzTol && maxIntensity > yThreshold then 
                        xData.Add currentMz
                        yData.Add currentIntensity              
                let finalX = xData.ToArray()
                let finalY = yData.ToArray()
                finalX, finalY                

                
        /// Returns a MzIntensityArray that containing the spectral centroids of the input spectra. 
        let toCentroid waveletF nScales yThreshold mzTol (*spacing_perc*) snrs_perc minSnr (refineMZ:bool) (sumIntensities:bool) (mzData: float []) (intensityData: float [])  =
            if mzData.Length < 3 then
                ([||], [||]) //, Array2D.empty
            else
            //FixedParameters
            let windowMaxLow = 5
            let windowMaxHigh = 5
            let initialWidthScaling = 1.; 
            let finalWidthScaling = 7.; 
            let incrementScaling = (finalWidthScaling-initialWidthScaling) / double(nScales-1);   
            let scalings = 
                Array.init nScales (fun index -> initialWidthScaling + float index * incrementScaling)
            // Calculation of the mzSpacing of each mz value
            let xSpacing = Care.createXspacing mzData 
            // Calculation of the average mz spacing of each mz value. This information is used to determine
            // the number of mz values for each individual mz value that are included in the computation of the individual
            // correlation
            let xSpacingAvg = Array.zeroCreate (mzData.Length) 
            let nPoints = Array3D.create 2 nScales mzData.Length 0
            getScales nScales scalings (*spacing_perc*) windowMaxLow windowMaxHigh mzData intensityData xSpacingAvg nPoints xSpacing    
            
            // Creates padded versions of the given data arrays
            let paddingPoints = 10000 // Can propably be reduced
            let mzDataPadded = Array.create (mzData.Length+2*paddingPoints) 0.
            let intensityDataPadded = Array.create (mzData.Length+2*paddingPoints) 0.   
            createPadding paddingPoints mzDataPadded intensityDataPadded mzData intensityData
            
            // Calculation of the wavelet correlation for each mzValue       
            let corrMatrixLength = 2*mzData.Length-1 //almost twice as long because also the correlation for Midpoints will be calculated.
            let corrMatrix = Array2D.create nScales corrMatrixLength 0. 
            calcCorrWaveletAndDataMatrix waveletF nScales paddingPoints mzDataPadded intensityDataPadded mzData intensityData nPoints xSpacingAvg scalings corrMatrix
            |> ignore 
            // Compares the previously selected correlation values of each mz value locally and selects the best scoring point
            let allLines = ResizeArray<RidgeLine>() 
            //
            getSNRFilteredPeakLines snrs_perc minSnr mzTol nScales mzData allLines corrMatrix

            refinePeaks refineMZ sumIntensities yThreshold mzTol scalings mzData intensityData allLines xSpacingAvg


        let toCentroidWithRicker2D (centroidParams:WaveletParameters) (mzData: float []) (intensityData: float [])  = 
            toCentroid ricker2d centroidParams.NumberOfScales centroidParams.YThreshold centroidParams.MzTolerance (*spacing_perc*) centroidParams.SNRS_Percentile centroidParams.MinSNR centroidParams.RefineMZ centroidParams.SumIntensities mzData intensityData


    module SecondDerivative =
        ///
        type SecondDerivativeParameters = {
            WindowSizeSGfilter    : int
            YThreshold            : float
            MzTolerance           : float
            SNRS_Windowsize       : int 
            SNRS_Percentile       : float
            RefineMZ              : bool
            SumIntensities        : bool
            }

        //toCentroid ricker2d nScales yYhreshold mzTol (*spacing_perc*) snrs_perc minSnr mzData intensityData
        ///
        let createSecondDerivativeParameters windowSizeSGfilter yThreshold mzTolerance sNRS_Windowsize sNRS_Percentile refineMZ sumIntensities = {
            WindowSizeSGfilter  = windowSizeSGfilter
            YThreshold          = yThreshold 
            MzTolerance         = mzTolerance
            SNRS_Windowsize     = sNRS_Windowsize
            SNRS_Percentile     = sNRS_Percentile
            RefineMZ            = refineMZ
            SumIntensities      = sumIntensities
            }

        let toCentroid refineMz sumIntensities windowSizeSGfilter yThreshold mzTol windowSizeSNR snrs_perc (mzData:float []) (intensityData:float []) =
            if mzData.Length < 5 || intensityData.Length < 5 then [||],[||]
            else
            // Step 1: Calculate negative snd derivative of the intensity Data
            let negSndDev = 
                FSharp.Stats.Signal.Filtering.savitzkyGolay windowSizeSGfilter 3 2 5 intensityData 
                |> Array.ofSeq
                |> Array.map (fun x -> x * -1.)

            // Step 2: Calculate noise threshold in snd derivative space;
            let minSnr = 1.
            let nNoiseBins = negSndDev.Length / windowSizeSNR + 1 //number of NoiseBins
            let noises = Array.create nNoiseBins 0.0

            for i = 0 to nNoiseBins-1 do
                let windowLow = 
                    match i * windowSizeSNR with 
                    | x when x = 0 -> 3
                    | x -> x  // inclusive
                let mutable windowHigh = 
                    let windowHigh' = windowLow + windowSizeSNR
                    if windowHigh' >= negSndDev.Length then 
                        negSndDev.Length
                    else
                        windowHigh'  
                if i = nNoiseBins-1 
                    then windowHigh <- negSndDev.Length
                       
                let unsortedData = ResizeArray<float>()

                for j = windowLow to windowHigh - 3 do
                     if (negSndDev.[j] < negSndDev.[j - 1] && negSndDev.[j] < negSndDev.[j + 1] 
                        && negSndDev.[j - 1] <= negSndDev.[j - 2] && negSndDev.[j + 1] <= negSndDev.[j + 2]) then unsortedData.Add (abs intensityData.[j])

                let nTot = unsortedData.Count
                let sortedData = unsortedData |> Seq.toArray |> Array.sort
                let mutable noise = Care.scoreAtPercentile snrs_perc nTot sortedData  //calcs Noisethreshold in bin i 
                if noise < 1.0 then noise <- 1.0
                noises.[i] <- noise
    
            // Step 3: find maxima above the noise level
            let filteredPeaks =
                let filteredPeaks = ResizeArray<(float*float)*int>()
                negSndDev
                |> (FSharp.Stats.Signal.PeakDetection.localMaximaIdx 0. mzData) 
                |> Array.iter (fun i -> 
                                let noiseBin = 
                                    let tmp = i/windowSizeSNR
                                    match tmp with
                                    | tmp when tmp > nNoiseBins-1 -> nNoiseBins-1 
                                    | _ -> tmp
                                let snr = intensityData.[i] / noises.[noiseBin]
                                if snr > minSnr then     
                                    filteredPeaks.Add ((mzData.[i],intensityData.[i]), i)
                               )
                filteredPeaks
                |> Array.ofSeq
            // Step 4: refine Peak positions to real maxima
            let refinedCentroids = 
                let xData = ResizeArray<float>()
                let yData = ResizeArray<float>()
                filteredPeaks
                |> Array.iteri (fun idx ((mz,intensity),i) -> 
                                    let startFittingPoint  = Care.getColLowBound mzData i (mzTol/2.)
                                    let endFittingPoint    = Care.getColHighBound mzData i (mzTol/2.)
                                    let mutable maxIntensityMZ = 0.0
                                    let mutable maxIntensity = 0.0 
                                    let mutable weightedMzSum = 0.0
                                    let mutable intensitySum = 0.0 
                                    for j = startFittingPoint to endFittingPoint do 
                                        intensitySum <- intensitySum + intensityData.[j] 
                                        weightedMzSum <- weightedMzSum + (mzData.[j] * intensityData.[j])
                                        if intensityData.[j] >= maxIntensity
                                            then
                                                maxIntensity   <- intensityData.[j]
                                                maxIntensityMZ <- mzData.[j] 
                                    let currentMz = 
                                        if refineMz then 
                                            let weightedMeanMz = weightedMzSum / intensitySum
                                            weightedMeanMz 
                                        else
                                            maxIntensityMZ
                                    let currentIntensity =
                                        if sumIntensities then 
                                            intensitySum
                                        else
                                            maxIntensity
                                    let xDataL = xData.Count    
                                    if xDataL = 0 && currentIntensity > yThreshold then 
                                        xData.Add (currentMz)
                                        yData.Add (currentIntensity)
                                    else
                                        let prevMz,prevMaxIntensity= xData.[xDataL-1],yData.[xDataL-1]
                                        if currentMz - prevMz >= mzTol && currentIntensity > yThreshold then
                                            xData.Add (currentMz)
                                            yData.Add (currentIntensity)
                                        elif currentMz  - prevMz <= mzTol && prevMaxIntensity < currentIntensity then
                                             xData.[xDataL-1] <- currentMz
                                             yData.[xDataL-1] <- currentIntensity
                               )
                xData.ToArray(),yData.ToArray()
            refinedCentroids

        let toCentroidSecondDerivative (centroidParams:SecondDerivativeParameters) (mzData: float []) (intensityData: float []) = 
            toCentroid centroidParams.RefineMZ centroidParams.SumIntensities centroidParams.WindowSizeSGfilter centroidParams.YThreshold centroidParams.MzTolerance centroidParams.SNRS_Windowsize centroidParams.SNRS_Percentile mzData intensityData
            
    /// Returns mzIntensityArray after noise reduction 
    let filterByIntensitySNR perc minSnr (mzData: float []) (intensityData:float [])  = 
        let noise = Care.scoreAtPercentile perc (intensityData.Length-1) intensityData 
        let filteredMz             = new ResizeArray<float>() 
        let filteredIntensityArray = new ResizeArray<float>()
        for i=0 to intensityData.Length-1 do
            let snr = intensityData.[i]/noise
            if snr > minSnr then
                filteredMz.Add mzData.[i]
                filteredIntensityArray.Add intensityData.[i]
        filteredMz |> Seq.toArray, filteredIntensityArray |> Seq.toArray

    /// Returns mzIntensityArray consisting of centroided Peaks. 
    let toCentroid centroidF (mzData:float[]) (intensityData:float[]) =
        if mzData.Length < 3  then 
            [||], [||]
        else
            centroidF mzData intensityData

    /// Returns mzIntensityArray consisting of centroided Peaks. 
    let windowToCentroid centroidF (mzData:float[]) (intensityData:float[]) lowerIdx upperIdx =
        if mzData.Length > 2  && lowerIdx>=0 && mzData.Length-2 >=upperIdx then 
            centroidF mzData.[lowerIdx.. upperIdx] intensityData.[lowerIdx.. upperIdx]  
        else [|float lowerIdx |], [|float upperIdx |]
        
    /// Returns a Index that accesses the mzData Array at a position determined by a precursorMz at the lower end of a given windowwidth 
    let lowerIdxBy (mzData:float []) windowWidth preCursorMZ =
        match Array.tryFindIndex (fun x -> x > preCursorMZ-(windowWidth/2.)) mzData with 
        | Some x -> x 
        | None -> mzData.Length-1 

    /// Returns a Index that accesses the mzData Array at a position determined by a precursorMz at the upper end of a given windowwidth 
    let upperIdxBy (mzData:float []) windowWidth preCursorMZ = 
        match Array.tryFindIndexBack (fun x -> x < preCursorMZ+(windowWidth/2.)) mzData with 
        | Some x -> x 
        | None  -> mzData.Length-1 
    
    // Returns mzIntensityArray consisting of centroided Peaks. 
    let windowToCentroidBy centroidF (mzData:float[]) (intensityData:float[]) windowWidth centerMass =
        if mzData.Length > 2 then
            let lowerIdx = lowerIdxBy mzData windowWidth centerMass
            let upperIdx = upperIdxBy mzData windowWidth centerMass
            windowToCentroid centroidF mzData intensityData lowerIdx upperIdx
        else [||], [||]



