System.h

Go to the documentation of this file.

#ifndef SimTK_SimTKCOMMON_SYSTEM_H_
#define SimTK_SimTKCOMMON_SYSTEM_H_

/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTKcommon                               *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK Core biosimulation toolkit originating from      *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2006-7 Stanford University and the Authors.         *
 * Authors: Michael Sherman                                                   *
 * Contributors: Peter Eastman                                                *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "SimTKcommon/basics.h"
#include "SimTKcommon/Simmatrix.h"
#include "SimTKcommon/internal/State.h"

namespace SimTK {

class Subsystem;
class DecorativeGeometry;
class DefaultSystemSubsystem;

SimTK_DEFINE_UNIQUE_INDEX_TYPE(EventId)


class SimTK_SimTKCOMMON_EXPORT System {
public:
    class Guts; // local; name is System::Guts
    friend class Guts;
private:
    // This is the only data member in this class. Also, any class derived from
    // System must have *NO* data members at all (data goes in the Guts class).
    Guts* guts;
public:
    System() : guts(0) { }
    System(const System&);
    System& operator=(const System&);
    ~System();

    const String& getName()    const;
    const String& getVersion() const;


    class ProjectOptions {
        unsigned long optionSet;
        explicit ProjectOptions(unsigned int o) : optionSet(o) { }
    public:

        enum Option {
            None   = 0x00,

            Q      = 0x01,
            U      = 0x02,
            QError = 0x04,
            UError = 0x08,

            PositionOnly = (Q|QError),
            VelocityOnly = (U|UError),
            All          = (PositionOnly|VelocityOnly)
        };

        ProjectOptions() : optionSet(0) { }

        // This is an implicit conversion
        ProjectOptions(Option opt) : optionSet((unsigned long)opt) { }

        // Implicit conversion to bool when needed
        operator bool() const {return optionSet != 0;}
        bool hasAnyPositionOptions() const {return (optionSet&(unsigned long)PositionOnly) != 0;}
        bool hasAnyVelocityOptions() const {return (optionSet&(unsigned long)VelocityOnly) != 0;}
        bool isEmpty() const {return optionSet==0;}

        bool isOptionSet(Option opt) const {return (optionSet&(unsigned long)opt) != 0;}
        void clear() {optionSet=0;}
        void clearOption(Option opt) {optionSet &= ~(unsigned long)opt;}
        void setOption  (Option opt) {optionSet |= (unsigned long)opt;}

        // Set operators: not, or, and, set difference
        ProjectOptions operator~() const {return ProjectOptions( (~optionSet) & (unsigned long)All );}
        ProjectOptions& operator|=(ProjectOptions opts) {optionSet |= opts.optionSet; return *this;}
        ProjectOptions& operator&=(ProjectOptions opts) {optionSet &= opts.optionSet; return *this;}
        ProjectOptions& operator-=(ProjectOptions opts) {optionSet &= ~opts.optionSet; return *this;}

        ProjectOptions& operator|=(Option opt) {setOption(opt); return *this;}
        ProjectOptions& operator-=(Option opt) {clearOption(opt); return *this;}
    };

        // STATISTICS //

    void resetAllCountersToZero();

        // Realization

    long getNumRealizationsOfThisStage(Stage) const;

    long getNumRealizeCalls() const;

        // Projection

    long getNumQProjections() const;
    long getNumUProjections() const;
    long getNumQErrorEstimateProjections() const;
    long getNumUErrorEstimateProjections() const;

    long getNumProjectCalls() const;

        // Event handling and reporting

    long getNumHandlerCallsThatChangedStage(Stage) const;

    long getNumHandleEventCalls() const;

    long getNumReportEventCalls() const;


        // REALIZATION //

    const State& realizeTopology() const;


    const State& getDefaultState() const;
    State&       updDefaultState();

    void realizeModel(State&) const;

    void realize(const State& s, Stage g = Stage::HighestValid) const;

    // Generate all decorative geometry computable at a specific stage. This will
    // throw an exception if the state hasn't already been realized
    // to that stage. Note that the list is not inclusive -- you have to
    // request geometry from each stage to get all of it.
    // This routine asks each subsystem in succession to generate its decorative geometry
    // and append it to the end of the vector.
    // If the stage is Stage::Topology, realizeTopology() must already have
    // been called but the State is ignored.
    void calcDecorativeGeometryAndAppend(const State&, Stage, std::vector<DecorativeGeometry>&) const;

        // THE CONTINUOUS SYSTEM //

            // UNCONSTRAINED

    Real calcTimescale(const State&) const;

    void calcYUnitWeights(const State&, Vector& weights) const;

            // CONSTRAINED

    void project(State&, Real consAccuracy, const Vector& yweights,
                 const Vector& ootols, Vector& yerrest, ProjectOptions=ProjectOptions::All) const;

    void calcYErrUnitTolerances(const State&, Vector& tolerances) const;


        // THE DISCRETE SYSTEM //

    class EventTriggerInfo;

    void setHasTimeAdvancedEvents(bool); // default=false
    bool hasTimeAdvancedEvents() const;


    enum EventCause {
        TriggeredEvents    =1,
        ScheduledEvents    =2,
        TimeAdvancedEvent  =3,
        TerminationEvent   =4,

        InvalidEventCause  = -1
    };
    static const char* getEventCauseName(EventCause);



    void handleEvents
       (State&, EventCause, const std::vector<EventId>& eventIds,
        Real accuracy, const Vector& yWeights, const Vector& ooConstraintTols,
        Stage& lowestModified, bool& shouldTerminate) const;
    
    
    void reportEvents(const State& s, EventCause cause, const std::vector<EventId>& eventIds) const;

    void calcEventTriggerInfo(const State&, std::vector<EventTriggerInfo>&) const;

    void calcTimeOfNextScheduledEvent(const State&, Real& tNextEvent, std::vector<EventId>& eventIds, bool includeCurrentTime) const;

    void calcTimeOfNextScheduledReport(const State&, Real& tNextEvent, std::vector<EventId>& eventIds, bool includeCurrentTime) const;
    
    //TODO: these operators should be provided by the Vector class where they
    //can be perfomed more efficiently.

    static Real calcWeightedRMSNorm(const Vector& values, const Vector& weights) {
        assert(weights.size() == values.size());
        if (values.size()==0) return 0;
        Real sumsq = 0;
        for (int i=0; i<values.size(); ++i) {
            const Real wv = weights[i]*values[i];
            sumsq += wv*wv;
        }
        return std::sqrt(sumsq/weights.size());
    }

    static Real calcWeightedInfinityNorm(const Vector& values, const Vector& weights) {
        assert(weights.size() == values.size());
        if (values.size()==0) return 0;
        Real maxval = 0;
        for (int i=0; i<values.size(); ++i) {
            const Real wv = std::abs(weights[i]*values[i]);
            if (wv > maxval) maxval=wv;
        }
        return maxval;
    }

    SubsystemIndex adoptSubsystem(Subsystem& child);

    int getNSubsystems() const;
    const Subsystem& getSubsystem(SubsystemIndex)   const;
    Subsystem&       updSubsystem(SubsystemIndex);
    const DefaultSystemSubsystem& getDefaultSubsystem() const;
    DefaultSystemSubsystem& updDefaultSubsystem();

    // Internal use only
    bool isOwnerHandle() const;
    bool isEmptyHandle() const;

    // There can be multiple handles on the same System.
    bool isSameSystem(const System& otherSystem) const;

    bool systemTopologyHasBeenRealized() const;

    // dynamic_cast the returned reference to a reference to your concrete Guts
    // class.
    const System::Guts& getSystemGuts() const {assert(guts); return *guts;}
    System::Guts&       updSystemGuts()       {assert(guts); return *guts;}

    // Put new *unowned* Guts into this *empty* handle and take over ownership.
    // If this handle is already in use, or if Guts is already owned this
    // routine will throw an exception.
    void adoptSystemGuts(System::Guts* g);

    explicit System(System::Guts* g) : guts(g) { }
    bool hasGuts() const {return guts!=0;}

private:
    class EventTriggerInfoRep;
};

inline static System::ProjectOptions operator|(System::ProjectOptions::Option  o1,    System::ProjectOptions::Option o2)    {return System::ProjectOptions(o1) |= o2;}
inline static System::ProjectOptions operator|(System::ProjectOptions          opts,  System::ProjectOptions::Option o)     {return opts |= o;}
inline static System::ProjectOptions operator|(System::ProjectOptions::Option  o,     System::ProjectOptions         opts)  {return opts |= o;}
inline static System::ProjectOptions operator&(System::ProjectOptions::Option  o1,    System::ProjectOptions::Option o2)    {return System::ProjectOptions(o1) &= o2;}
inline static System::ProjectOptions operator&(System::ProjectOptions          opts,  System::ProjectOptions::Option o)     {return opts &= o;}
inline static System::ProjectOptions operator&(System::ProjectOptions::Option  o,     System::ProjectOptions         opts)  {return opts &= o;}
inline static System::ProjectOptions operator~(System::ProjectOptions::Option  o)                                           {return ~System::ProjectOptions(o);}
inline static System::ProjectOptions operator-(System::ProjectOptions          opts,  System::ProjectOptions::Option o)     {return opts -= o;}
inline static System::ProjectOptions operator-(System::ProjectOptions          opts1, System::ProjectOptions         opts2) {return opts1 -= opts2;}



class SimTK_SimTKCOMMON_EXPORT System::EventTriggerInfo {
public:
    EventTriggerInfo();
    explicit EventTriggerInfo(EventId eventId);
    ~EventTriggerInfo();
    EventTriggerInfo(const EventTriggerInfo&);
    EventTriggerInfo& operator=(const EventTriggerInfo&);

    EventId getEventId() const; // returns -1 if not set
    bool shouldTriggerOnRisingSignTransition()  const; // default=true
    bool shouldTriggerOnFallingSignTransition() const; // default=true
    Real getRequiredLocalizationTimeWindow()    const; // default=0.1

    // These return the modified 'this', like assignment operators.
    EventTriggerInfo& setEventId(EventId);
    EventTriggerInfo& setTriggerOnRisingSignTransition(bool);
    EventTriggerInfo& setTriggerOnFallingSignTransition(bool);
    EventTriggerInfo& setRequiredLocalizationTimeWindow(Real);

    // TODO: switch to SignTransitionSet and trash EventStatus altogether.
    EventStatus::EventTrigger calcTransitionMask() const {
        unsigned mask = 0;
        if (shouldTriggerOnRisingSignTransition()) {
            mask |= EventStatus::NegativeToPositive;
        }
        if (shouldTriggerOnFallingSignTransition()) {
            mask |= EventStatus::PositiveToNegative;
        }
        return EventStatus::EventTrigger(mask);
    }

    EventStatus::EventTrigger calcTransitionToReport
       (EventStatus::EventTrigger transitionSeen) const
    {
        // report -1 to 1 or 1 to -1 as appropriate
        if (transitionSeen & EventStatus::Rising)
            return EventStatus::NegativeToPositive;
        if (transitionSeen & EventStatus::Falling)
            return EventStatus::PositiveToNegative;
        assert(!"impossible event transition situation");
        return EventStatus::NoEventTrigger;
    }

private:
    // opaque implementation for future binary compatibility
    System::EventTriggerInfoRep* rep;

    const System::EventTriggerInfoRep& getRep() const {assert(rep); return *rep;}
    System::EventTriggerInfoRep&       updRep()       {assert(rep); return *rep;}
};

} // namespace SimTK

#endif // SimTK_SimTKCOMMON_SYSTEM_H_