TmMolStar/src/mol-geo/representation/structure/visual/util/polymer.ts

/**
 * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
 *
 * @author Alexander Rose <alexander.rose@weirdbyte.de>
 */

import { Unit, Element, StructureProperties, Model } from 'mol-model/structure';
import { Segmentation, OrderedSet, Interval } from 'mol-data/int';
import { MoleculeType, SecondaryStructureType } from 'mol-model/structure/model/types';
import Iterator from 'mol-data/iterator';
import { Vec3 } from 'mol-math/linear-algebra';
import { SymmetryOperator } from 'mol-math/geometry';
import SortedRanges from 'mol-data/int/sorted-ranges';

export function getPolymerRanges(unit: Unit) {
    switch (unit.kind) {
        case Unit.Kind.Atomic: return unit.model.atomicHierarchy.polymerRanges
        case Unit.Kind.Spheres: return unit.model.coarseHierarchy.spheres.polymerRanges
        case Unit.Kind.Gaussians: return unit.model.coarseHierarchy.gaussians.polymerRanges
    }
}

export function getPolymerElementCount(unit: Unit) {
    let count = 0
    const { elements } = unit
    const polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), elements)
    switch (unit.kind) {
        case Unit.Kind.Atomic:
            const residueIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueSegments, elements)
            while (polymerIt.hasNext) {
                const polymerSegment = polymerIt.move()
                residueIt.setSegment(polymerSegment)
                while (residueIt.hasNext) {
                    const residueSegment = residueIt.move()
                    const { start, end } = residueSegment
                    if (OrderedSet.areIntersecting(Interval.ofBounds(elements[start], elements[end - 1]), elements)) ++count
                }
            }
            break
        case Unit.Kind.Spheres:
        case Unit.Kind.Gaussians:
            while (polymerIt.hasNext) {
                const { start, end } = polymerIt.move()
                // TODO add OrderedSet.intersectionSize
                count += OrderedSet.size(OrderedSet.intersect(Interval.ofBounds(elements[start], elements[end - 1]), elements))
            }
            break
    }
    return count
}

function getTraceName(l: Element.Location) {
    const compId = StructureProperties.residue.label_comp_id(l)
    const chemCompMap = l.unit.model.properties.chemicalComponentMap
    const cc = chemCompMap.get(compId)
    const moleculeType = cc ? cc.moleculeType : MoleculeType.unknown
    let traceName = ''
    if (moleculeType === MoleculeType.protein) {
        traceName = 'CA'
    } else if (moleculeType === MoleculeType.DNA || moleculeType === MoleculeType.RNA) {
        traceName = 'P'
    }
    return traceName
}

function setTraceElement(l: Element.Location, residueSegment: Segmentation.Segment<Element>) {
    const elements = l.unit.elements
    l.element = elements[residueSegment.start]
    const traceName = getTraceName(l)

    for (let j = residueSegment.start, _j = residueSegment.end; j < _j; j++) {
        l.element = elements[j];
        if (StructureProperties.atom.label_atom_id(l) === traceName) return j
    }
    return residueSegment.end - 1
}


function getTraceName2(model: Model, residueModelIndex: number) {
    const compId = model.atomicHierarchy.residues.label_comp_id.value(residueModelIndex)
    const chemCompMap = model.properties.chemicalComponentMap
    const cc = chemCompMap.get(compId)
    const moleculeType = cc ? cc.moleculeType : MoleculeType.unknown
    let traceName = ''
    if (moleculeType === MoleculeType.protein) {
        traceName = 'CA'
    } else if (moleculeType === MoleculeType.DNA) {
        // traceName = 'P'
        traceName = 'C3\''
    } else if (moleculeType === MoleculeType.RNA) {
        // traceName = 'P'
        traceName = 'C4\''
    }
    return traceName
}

function getTraceElement2(model: Model, residueModelSegment: Segmentation.Segment<Element>) {
    const traceName = getTraceName2(model, residueModelSegment.index)

    for (let j = residueModelSegment.start, _j = residueModelSegment.end; j < _j; j++) {
        if (model.atomicHierarchy.atoms.label_atom_id.value(j) === traceName) return j
    }
    console.log(`trace name element "${traceName}" not found`, { ...residueModelSegment })
    return residueModelSegment.start
}

function getDirectionName2(model: Model, residueModelIndex: number) {
    const compId = model.atomicHierarchy.residues.label_comp_id.value(residueModelIndex)
    const chemCompMap = model.properties.chemicalComponentMap
    const cc = chemCompMap.get(compId)
    const moleculeType = cc ? cc.moleculeType : MoleculeType.unknown
    let traceName = ''
    if (moleculeType === MoleculeType.protein) {
        traceName = 'O'
    } else if (moleculeType === MoleculeType.DNA) {
        traceName = 'O4\''
    } else if (moleculeType === MoleculeType.RNA) {
        traceName = 'C3\''
    } else {
        console.log('molecule type unknown', Number.isInteger(residueModelIndex) ? compId : residueModelIndex, chemCompMap)
    }
    return traceName
}

// function residueLabel(model: Model, rI: number) {
//     const { residues, chains, residueSegments, chainSegments } = model.atomicHierarchy
//     const { label_comp_id, label_seq_id } = residues
//     const { label_asym_id } = chains
//     const cI = chainSegments.segmentMap[residueSegments.segments[rI]]
//     return `${label_asym_id.value(cI)} ${label_comp_id.value(rI)} ${label_seq_id.value(rI)}`
// }

function getDirectionElement2(model: Model, residueModelSegment: Segmentation.Segment<Element>) {
    const traceName = getDirectionName2(model, residueModelSegment.index)

    for (let j = residueModelSegment.start, _j = residueModelSegment.end; j < _j; j++) {
        if (model.atomicHierarchy.atoms.label_atom_id.value(j) === traceName) return j
    }
    // console.log('direction name element not found', { ...residueModelSegment })
    return residueModelSegment.start
}


/** Iterates over consecutive pairs of residues/coarse elements in polymers */
export function PolymerBackboneIterator(unit: Unit): Iterator<PolymerBackbonePair> {
    switch (unit.kind) {
        case Unit.Kind.Atomic: return new AtomicPolymerBackboneIterator(unit)
        case Unit.Kind.Spheres:
        case Unit.Kind.Gaussians:
            return new CoarsePolymerBackboneIterator(unit)
    }
}

interface PolymerBackbonePair {
    centerA: Element.Location
    centerB: Element.Location
    indexA: number
    indexB: number
    posA: Vec3
    posB: Vec3
}

function createPolymerBackbonePair (unit: Unit) {
    return {
        centerA: Element.Location(unit),
        centerB: Element.Location(unit),
        indexA: 0,
        indexB: 0,
        posA: Vec3.zero(),
        posB: Vec3.zero()
    }
}

const enum AtomicPolymerBackboneIteratorState { nextPolymer, firstResidue, nextResidue }

export class AtomicPolymerBackboneIterator<T extends number = number> implements Iterator<PolymerBackbonePair> {
    private value: PolymerBackbonePair

    private polymerIt: SortedRanges.Iterator<Element>
    private residueIt: Segmentation.Iterator<Element>
    private polymerSegment: Segmentation.Segment<Element>
    private state: AtomicPolymerBackboneIteratorState = AtomicPolymerBackboneIteratorState.nextPolymer
    private pos: SymmetryOperator.CoordinateMapper

    hasNext: boolean = false;

    move() {
        const { residueIt, polymerIt, value, pos } = this

        if (this.state === AtomicPolymerBackboneIteratorState.nextPolymer) {
            while (polymerIt.hasNext) {
                this.polymerSegment = polymerIt.move();
                // console.log('polymerSegment', this.polymerSegment)
                residueIt.setSegment(this.polymerSegment);

                const residueSegment = residueIt.move();
                // console.log('first residueSegment', residueSegment, residueIt.hasNext)
                if (residueIt.hasNext) {
                    value.indexB = setTraceElement(value.centerB, residueSegment)
                    pos(value.centerB.element, value.posB)
                    this.state = AtomicPolymerBackboneIteratorState.nextResidue
                    break
                }
            }
        }

        if (this.state === AtomicPolymerBackboneIteratorState.nextResidue) {
            const residueSegment = residueIt.move();
            // console.log('next residueSegment', residueSegment)
            value.centerA.element = value.centerB.element
            value.indexA = value.indexB
            Vec3.copy(value.posA, value.posB)
            value.indexB = setTraceElement(value.centerB, residueSegment)
            pos(value.centerB.element, value.posB)

            if (!residueIt.hasNext) {
                // TODO need to advance to a polymer that has two or more residues (can't assume it has)
                this.state = AtomicPolymerBackboneIteratorState.nextPolymer
            }
        }

        this.hasNext = residueIt.hasNext || polymerIt.hasNext
        
        // console.log('hasNext', this.hasNext)
        // console.log('value', this.value)

        return this.value;
    }

    constructor(unit: Unit.Atomic) {
        const { residueSegments } = unit.model.atomicHierarchy
        // console.log('unit.elements', OrderedSet.toArray(unit.elements))
        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements)
        this.residueIt = Segmentation.transientSegments(residueSegments, unit.elements)
        this.pos = unit.conformation.invariantPosition
        this.value = createPolymerBackbonePair(unit)

        this.hasNext = this.residueIt.hasNext && this.polymerIt.hasNext
    }
}

const enum CoarsePolymerBackboneIteratorState { nextPolymer, nextElement }

export class CoarsePolymerBackboneIterator<T extends number = number> implements Iterator<PolymerBackbonePair> {
    private value: PolymerBackbonePair

    private polymerIt: SortedRanges.Iterator<Element>
    private polymerSegment: Segmentation.Segment<Element>
    private state: CoarsePolymerBackboneIteratorState = CoarsePolymerBackboneIteratorState.nextPolymer
    private pos: SymmetryOperator.CoordinateMapper
    private elementIndex: number

    hasNext: boolean = false;

    move() {
        const { polymerIt, value, pos } = this

        if (this.state === CoarsePolymerBackboneIteratorState.nextPolymer) {
            if (polymerIt.hasNext) {
                this.polymerSegment = polymerIt.move();
                console.log('polymer', this.polymerSegment)
                this.elementIndex = this.polymerSegment.start
                this.elementIndex += 1
                if (this.elementIndex + 1 < this.polymerSegment.end) {
                    value.centerB.element = value.centerB.unit.elements[this.elementIndex]
                    value.indexB = this.elementIndex
                    pos(value.centerB.element, value.posB)

                    this.state = CoarsePolymerBackboneIteratorState.nextElement
                } else {
                    this.state = CoarsePolymerBackboneIteratorState.nextPolymer
                }
            }
        }

        if (this.state === CoarsePolymerBackboneIteratorState.nextElement) {
            this.elementIndex += 1
            value.centerA.element = value.centerB.element
            value.indexA = value.indexB
            Vec3.copy(value.posA, value.posB)

            value.centerB.element = value.centerB.unit.elements[this.elementIndex]
            value.indexB = this.elementIndex
            pos(value.centerB.element, value.posB)

            if (this.elementIndex + 1 >= this.polymerSegment.end) {
                this.state = CoarsePolymerBackboneIteratorState.nextPolymer
            }
        }

        this.hasNext = this.elementIndex + 1 < this.polymerSegment.end || polymerIt.hasNext
        return this.value;
    }

    constructor(unit: Unit.Spheres | Unit.Gaussians) {
        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements);

        this.pos = unit.conformation.invariantPosition
        this.value = createPolymerBackbonePair(unit)

        console.log('CoarsePolymerBackboneIterator', this.polymerIt.hasNext)

        this.hasNext = this.polymerIt.hasNext
    }
}




/**
 * Iterates over individual residues/coarse elements in polymers of a unit while
 * providing information about the neighbourhood in the underlying model for drawing splines
 */
export function PolymerTraceIterator(unit: Unit): Iterator<PolymerTraceElement> {
    switch (unit.kind) {
        case Unit.Kind.Atomic: return new AtomicPolymerTraceIterator(unit)
        case Unit.Kind.Spheres:
        case Unit.Kind.Gaussians:
            return new CoarsePolymerTraceIterator(unit)
    }
}

interface PolymerTraceElement {
    center: Element.Location
    index: number
    first: boolean
    last: boolean
    secStrucType: SecondaryStructureType

    t0: Vec3
    t1: Vec3
    t2: Vec3
    t3: Vec3
    t4: Vec3

    d12: Vec3
    d23: Vec3
}

function createPolymerTraceElement (unit: Unit): PolymerTraceElement {
    return {
        center: Element.Location(unit),
        index: 0,
        first: false,
        last: false,
        secStrucType: SecondaryStructureType.create(SecondaryStructureType.Flag.NA),

        t0: Vec3.zero(),
        t1: Vec3.zero(),
        t2: Vec3.zero(),
        t3: Vec3.zero(),
        t4: Vec3.zero(),

        d12: Vec3.zero(),
        d23: Vec3.zero(),
    }
}

const enum AtomicPolymerTraceIteratorState { nextPolymer, nextResidue }

function setSegment (outSegment: Segmentation.Segment<Element>, index: number, segments: Segmentation<Element>, min: number, max: number): Segmentation.Segment<Element> {
    // index = Math.min(Math.max(0, index), segments.segments.length - 2)
    const _index = Math.min(Math.max(min, index), max)
    if (isNaN(_index)) console.log(_index, index, min, max)
    outSegment.index = _index
    outSegment.start = segments.segments[_index]
    outSegment.end = segments.segments[_index + 1]
    // console.log(index, {...outSegment}, {...boundingSegment}, segments.segments[boundingSegment.index])
    return outSegment
}

export class AtomicPolymerTraceIterator<T extends number = number> implements Iterator<PolymerTraceElement> {
    private value: PolymerTraceElement

    private polymerIt: SortedRanges.Iterator<Element>
    private residueIt: Segmentation.Iterator<Element>
    private residueSegmentMin: number
    private residueSegmentMax: number
    private state: AtomicPolymerTraceIteratorState = AtomicPolymerTraceIteratorState.nextPolymer
    private residueSegments: Segmentation<Element>

    private tmpSegment: Segmentation.Segment<Element>

    private unit: Unit.Atomic

    hasNext: boolean = false;

    private pos(target: Vec3, index: number) {
        target[0] = this.unit.model.atomicConformation.x[index]
        target[1] = this.unit.model.atomicConformation.y[index]
        target[2] = this.unit.model.atomicConformation.z[index]
    }

    updateResidueSegmentRange(polymerSegment: Segmentation.Segment<Element>) {
        const { polymerRanges, residueSegments } = this.unit.model.atomicHierarchy
        const sMin = polymerRanges[polymerSegment.index * 2]
        const sMax = polymerRanges[polymerSegment.index * 2 + 1]
        this.residueSegmentMin = residueSegments.segmentMap[sMin]
        this.residueSegmentMax = residueSegments.segmentMap[sMax]
    }

    move() {
        const { residueIt, polymerIt, value } = this

        if (this.state === AtomicPolymerTraceIteratorState.nextPolymer) {
            while (polymerIt.hasNext) {
                const polymerSegment = polymerIt.move();
                // console.log('polymerSegment', {...polymerSegment})
                residueIt.setSegment(polymerSegment);
                this.updateResidueSegmentRange(polymerSegment)
                if (residueIt.hasNext) {
                    this.state = AtomicPolymerTraceIteratorState.nextResidue
                    break
                }
            }
        }

        if (this.state === AtomicPolymerTraceIteratorState.nextResidue) {
            const { tmpSegment, residueSegments, residueSegmentMin, residueSegmentMax } = this
            const residueSegment = residueIt.move();
            const resSegIdx = residueSegment.index
            // console.log(residueLabel(this.unit.model, resSegIdx), resSegIdx, this.unit.model.properties.secondaryStructure.type[resSegIdx])
            value.index = setTraceElement(value.center, residueSegment)

            setSegment(tmpSegment, resSegIdx - 2, residueSegments, residueSegmentMin, residueSegmentMax)
            this.pos(value.t0, getTraceElement2(this.unit.model, tmpSegment))

            setSegment(tmpSegment, resSegIdx - 1, residueSegments, residueSegmentMin, residueSegmentMax)
            this.pos(value.t1, getTraceElement2(this.unit.model, tmpSegment))
            this.pos(value.d12, getDirectionElement2(this.unit.model, tmpSegment))

            setSegment(tmpSegment, resSegIdx, residueSegments, residueSegmentMin, residueSegmentMax)
            value.secStrucType = this.unit.model.properties.secondaryStructure.type[resSegIdx]
            this.pos(value.t2, getTraceElement2(this.unit.model, tmpSegment))
            this.pos(value.d23, getDirectionElement2(this.unit.model, tmpSegment))

            setSegment(tmpSegment, resSegIdx + 1, residueSegments, residueSegmentMin, residueSegmentMax)
            this.pos(value.t3, getTraceElement2(this.unit.model, tmpSegment))

            setSegment(tmpSegment, resSegIdx + 2, residueSegments, residueSegmentMin, residueSegmentMax)
            this.pos(value.t4, getTraceElement2(this.unit.model, tmpSegment))

            value.first = resSegIdx === residueSegmentMin
            value.last = resSegIdx === residueSegmentMax

            if (!residueIt.hasNext) {
                this.state = AtomicPolymerTraceIteratorState.nextPolymer
            }
        }

        this.hasNext = residueIt.hasNext || polymerIt.hasNext

        return this.value;
    }

    constructor(unit: Unit.Atomic) {
        const { residueSegments } = unit.model.atomicHierarchy
        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements)
        this.residueIt = Segmentation.transientSegments(residueSegments, unit.elements);
        this.residueSegments = residueSegments
        this.value = createPolymerTraceElement(unit)
        this.hasNext = this.residueIt.hasNext && this.polymerIt.hasNext

        this.tmpSegment = { index: 0, start: 0 as Element, end: 0 as Element }

        this.unit = unit
    }
}

export class CoarsePolymerTraceIterator<T extends number = number> implements Iterator<PolymerTraceElement> {
    private value: PolymerTraceElement

    hasNext: boolean = false;

    move() {
        return this.value;
    }

    constructor(unit: Unit.Spheres | Unit.Gaussians) {
        this.value = createPolymerTraceElement(unit)
        this.hasNext = false
    }
}