Subversion Repositories public iLand

/********************************************************************************************

**    iLand - an individual based forest landscape and disturbance model

**    http://iland.boku.ac.at

**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl

**

**    This program is free software: you can redistribute it and/or modify

**    it under the terms of the GNU General Public License as published by

**    the Free Software Foundation, either version 3 of the License, or

**    (at your option) any later version.

**

**    This program is distributed in the hope that it will be useful,

**    but WITHOUT ANY WARRANTY; without even the implied warranty of

**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

**    GNU General Public License for more details.

**

**    You should have received a copy of the GNU General Public License

**    along with this program.  If not, see <http://www.gnu.org/licenses/>.

********************************************************************************************/

#include "abe_global.h"

#include "fmstand.h"

#include "fmunit.h"

#include "management.h"

#include "fmtreelist.h"

#include "forestmanagementengine.h"

#include "mapgrid.h"

#include "fmstp.h"

#include "scheduler.h"

#include "fomescript.h"

#include "agent.h"

#include "agenttype.h"

#include "tree.h"

#include "species.h"

#include "statdata.h"

#include "debugtimer.h"

namespace ABE {

/** @class FMStand

    @ingroup abe

    The FMStand class encapsulates forest stands which are defined as polygons. FMStand tracks properties of the stands (e.g. mean volume), and

    is a central player in the ABE system.

  */

FMStand::FMStand(FMUnit *unit, const int id)

{

    mUnit = unit;

    mId = id;

    mInitialId = id;

    mPhase = Activity::Invalid;

    // testing:

    mPhase = Activity::Tending;

    mStandType = 1; // just testing...

    mU = 0, mSpeciesCompositionIndex = -1, mThinningIntensityClass = -1;

    newRotatation();

    mSTP = 0;

    mVolume = 0.;

    mAge = 0.;

    mTotalBasalArea = 0.;

    mStems = 0.;

    mDbh = 0.;

    mHeight = 0.;

    mScheduledHarvest = 0.;

    mFinalHarvested = 0.;

    mThinningHarvest = 0.;

    mDisturbed = 0.;

    mRotationStartYear = 0;

    mLastUpdate = -1.;

    mLastExecution = -1.;

    mCurrentIndex=-1;

    mLastExecutedIndex=-1;

    mLastRotationAge = -1;

    mArea = ForestManagementEngine::standGrid()->area(mId)/cRUArea;

}

void FMStand::initialize()

{

    if (!mSTP)

        throw IException(QString("FMStand::initialize, no valid STP for stand %1").arg(id()));

    // copy activity flags

    mStandFlags = mSTP->defaultFlags();

    mCurrentIndex=-1;

    mLastExecutedIndex=-1;

    mYearsToWait=0;

    mContextStr = QString("S%2Y%1:").arg(ForestManagementEngine::instance()->currentYear()).arg(id()); // initialize...

    // load data and aggregate averages

    reload();

    if (mRotationStartYear==0.) // only set if not explicitely set previously.

        mRotationStartYear = ForestManagementEngine::instance()->currentYear() - age();

    // when a stand is initialized, we assume that 20% of the standing volume

    // have been removed already.

    mRemovedVolumeTotal = volume() * 0.2;

    if (absoluteAge()>0)

        mMAItotal = volume() * 1.2 / absoluteAge();

    else

        mMAItotal = 0.;

    mMAIdecade = mMAItotal;

    mLastMAIVolume = volume();

    // find out the first activity...

    int min_years_to_wait = 100000;

    for (int i=0;i<mStandFlags.count(); ++i) {

        // run the onSetup event

        // specifically set 'i' as the activity to be evaluated.

        FomeScript::setExecutionContext(this);

        FomeScript::bridge()->activityObj()->setActivityIndex(i);

        mStandFlags[i].activity()->events().run(QStringLiteral("onSetup"), 0);

        if (!mStandFlags[i].enabled() || !mStandFlags[i].active())

            continue;

        // set active to false which have already passed

        if (!mStandFlags[i].activity()->isRepeatingActivity()) {

            if (!mStandFlags[i].activity()->schedule().absolute && mStandFlags[i].activity()->latestSchedule(U()) < age()) {

                mStandFlags[i].setActive(false);

            } else {

                int delta = mStandFlags[i].activity()->earlistSchedule(U()) - age();

                if (mStandFlags[i].activity()->schedule().absolute)

                    delta += age(); // absolute timing: starting from 0

                if (delta<min_years_to_wait) {

                    min_years_to_wait = qMax(delta,0); // limit to 0 years

                    mCurrentIndex = i; // first activity to execute

                }

            }

        }

    }

    if (mCurrentIndex==-1) {

        // the stand is "outside" the time frames provided by the activities.

        // set the last activity with "force" = true as the active

        for (int i=mStandFlags.count()-1;i>=0; --i)

            if (mStandFlags[i].enabled() && mStandFlags[i].activity()->schedule().force_execution==true) {

                mCurrentIndex = i;

                break;

            }

    }

    if (min_years_to_wait<100000)

        sleep(min_years_to_wait);

    // call onInit handler on the level of the STP

    mSTP->events().run(QStringLiteral("onInit"), this);

    if (mCurrentIndex>-1) {

        mStandFlags[mCurrentIndex].activity()->events().run(QStringLiteral("onEnter"), this);

        // if it is a scheduled activity, then execute (to get initial estimates for harvests)

        if (currentFlags().isScheduled())

            executeActivity(currentActivity());

    }

}

void FMStand::reset(FMSTP *stp)

{

    mSTP = stp;

    newRotatation();

    mCurrentIndex = -1;

}

void FMStand::checkArea()

{

    mArea = ForestManagementEngine::standGrid()->area(mId)/cRUArea;

}

bool relBasalAreaIsHigher(const SSpeciesStand &a, const SSpeciesStand &b)

{

    return a.relBasalArea > b.relBasalArea;

}

void FMStand::reload(bool force)

{

    if (!force && mLastUpdate == ForestManagementEngine::instance()->currentYear())

        return;

    DebugTimer t("ABE:FMStand::reload");

    // load all trees that are located on this stand

    mTotalBasalArea = 0.;

    mVolume = 0.;

    mAge = 0.;

    mStems = 0.;

    mDbh = 0.;

    mHeight = 0.;

    mTopHeight = 0.;

    mLastUpdate = ForestManagementEngine::instance()->currentYear();

    mSpeciesData.clear();

    // load all trees of the forest stand (use the treelist of the current execution context)

    FMTreeList *trees = ForestManagementEngine::instance()->scriptBridge()->treesObj();

    trees->setStand(this);

    trees->loadAll();

    //qDebug() << "fmstand-reload: load trees from map:" << t.elapsed();

    // use: value_per_ha = value_stand * area_factor

    double area_factor = 1. / area();

    const QVector<QPair<Tree*, double> > &treelist = trees->trees();

    // calculate top-height: diameter of the 100 thickest trees per ha

    QVector<double> dbhvalues;

    dbhvalues.reserve(trees->trees().size());

    for ( QVector<QPair<Tree*, double> >::const_iterator it=treelist.constBegin(); it!=treelist.constEnd(); ++it)

        dbhvalues.push_back(it->first->dbh());

    double topheight_threshhold=0.;

    double topheight_height = 0.;

    int topheight_trees = 0;

    if (treelist.size()>0) {

        StatData s(dbhvalues);

        topheight_threshhold= s.percentile( 100*(1- area()*100/treelist.size()) ); // sorted ascending -> thick trees at the end of the list

    }

    for ( QVector<QPair<Tree*, double> >::const_iterator it=treelist.constBegin(); it!=treelist.constEnd(); ++it) {

        double ba = it->first->basalArea() * area_factor;

        mTotalBasalArea+=ba;

        mVolume += it->first->volume() * area_factor;

        mAge += it->first->age()*ba;

        mDbh += it->first->dbh()*ba;

        mHeight += it->first->height()*ba;

        mStems++;

        SSpeciesStand &sd = speciesData(it->first->species());

        sd.basalArea += ba;

        if (it->first->dbh() >= topheight_threshhold) {

            topheight_height += it->first->height();

            ++topheight_trees;

        }

    }

    if (mTotalBasalArea>0.) {

        mAge /= mTotalBasalArea;

        mDbh /= mTotalBasalArea;

        mHeight /= mTotalBasalArea;

        for (int i=0;i<mSpeciesData.count();++i) {

            mSpeciesData[i].relBasalArea =  mSpeciesData[i].basalArea / mTotalBasalArea;

        }

    }

    if (topheight_trees>0) {

        mTopHeight = topheight_height / double(topheight_trees);

    }

    mStems *= area_factor; // convert to stems/ha

    // sort species data by relative share....

    std::sort(mSpeciesData.begin(), mSpeciesData.end(), relBasalAreaIsHigher);

}

// return stand area in ha

double FMStand::absoluteAge() const

{

    return ForestManagementEngine::instance()->currentYear() - mRotationStartYear;

}

bool FMStand::execute()

{

    //  the age of the stand increases by one

    mAge++;

    // do nothing if we are still waiting (sleep)

    if (mYearsToWait>0) {

        if (--mYearsToWait > 0) {

            return false;

        }

    }

    mContextStr = QString("S%2Y%1:").arg(ForestManagementEngine::instance()->currentYear()).arg(id());

    // what to do if there is no active activity??

    if (mCurrentIndex==-1) {

        if (trace())

            qCDebug(abe) << context() << "*** No action - no currently active activity ***";

        return false;

    }

    if (trace())

        qCDebug(abe) << context() << "*** start evaulate activity:" << currentActivity()->name();

    // do nothing if there is already an activity in the scheduler

    if (currentFlags().isPending()) {

        if (trace())

            qCDebug(abe) << context() << "*** No action - stand in the scheduler. ***";

        return false;

    }

    // do nothing if the the current year is not within the window of opportunity of the activity

    double p_schedule = currentActivity()->scheduleProbability(this);

    if (p_schedule == -1.) {

        if (trace())

            qCDebug(abe)<< context()  << "*** Activity expired. ***";

        // cancel the activity

        currentFlags().setActive(false);

        afterExecution(true);

        return false;

    }

    if (p_schedule>=0. && p_schedule < 0.00001) {

        if (trace())

            qCDebug(abe)<< context()  << "*** No action - Schedule probability 0. ***";

        return false;

    }

    // we need to renew the stand data

    reload();

    // check if there are some constraints that prevent execution....

    double p_execute = currentActivity()->execeuteProbability(this);

    if (p_execute == 0.) {

        if (trace())

            qCDebug(abe)<< context() << "*** No action - Constraints preventing execution. ***";

        return false;

    }

    // ok, we should execute the current activity.

    // if it is not scheduled, it is executed immediately, otherwise a ticket is created.

    if (currentFlags().isScheduled()) {

        // ok, we schedule the current activity

        if (trace())

            qCDebug(abe)<< context() << "adding ticket for execution.";

        mScheduledHarvest = 0.;

        bool should_schedule = currentActivity()->evaluate(this);

        if (trace())

            qCDebug(abe) << context() << "evaluated stand. add a ticket:" << should_schedule;

        if (should_schedule) {

            mUnit->scheduler()->addTicket(this, &currentFlags(), p_schedule, p_execute );

        } else {

            // cancel the activity

            currentFlags().setActive(false);

            afterExecution(true);

        }

        return should_schedule;

    } else {

        // execute immediately

        if (trace())

            qCDebug(abe) << context() << "executing activty" << currentActivity()->name();

        mScheduledHarvest = 0.;

        bool executed = currentActivity()->execute(this);

        if (!currentActivity()) // special case: the activity invalidated the active activtity

            return executed;

        if (!currentActivity()->isRepeatingActivity()) {

            currentFlags().setActive(false);

            afterExecution(!executed); // check what comes next for the stand

        }

        return executed;

    }

}

bool FMStand::executeActivity(Activity *act)

{

    int old_activity_index = mCurrentIndex;

    int new_index = stp()->activityIndex(act);

    bool result=false;

    if (new_index>-1) {

        mCurrentIndex = new_index;

        int old_years = mYearsToWait;

        mYearsToWait = 0;

        result = execute();

        mAge--; // undo modification of age

        mYearsToWait = old_years; // undo...

    }

    mCurrentIndex = old_activity_index;

    return result;

}

bool FMStand::afterExecution(bool cancel)

{

    // check if an agent update is necessary

    unit()->agent()->type()->agentUpdateForStand(this, currentFlags().activity()->name(), -1);

    // is called after an activity has run

    int tmin = 10000000;

    int indexmin = -1;

    for (int i=0;i<mStandFlags.count(); ++i) {

        if (mStandFlags[i].isForcedNext()) {

            mStandFlags[i].setForceNext(false); // reset flag

            indexmin = i;

            break; // we "jump" to this activity

        }

    }

    if (indexmin == -1) {

        // check if a restart is needed

        // TODO: find a better way!!

        if (currentFlags().isFinalHarvest()) {

            // we have reached the last activity

            for (int i=0;i<mStandFlags.count(); ++i)

                mStandFlags[i].setActive(true);

            newRotatation();

            reload();

        }

        // look for the next (enabled) activity.

        for (int i=0;i<mStandFlags.count(); ++i) {

            if ( mStandFlags[i].enabled() && mStandFlags[i].active() && !mStandFlags[i].isRepeating())

                if (mStandFlags[i].activity()->earlistSchedule() < tmin) {

                    tmin =  mStandFlags[i].activity()->earlistSchedule();

                    indexmin = i;

                }

        }

    }

    if (!cancel)

        currentActivity()->events().run(QStringLiteral("onExecuted"),this);

    else

        currentActivity()->events().run(QStringLiteral("onCancel"),this);

    if (indexmin != mCurrentIndex) {

        // call events:

        currentActivity()->events().run(QStringLiteral("onExit"), this);

        if (indexmin>-1 && indexmin<mStandFlags.count())

            mStandFlags[indexmin].activity()->events().run(QStringLiteral("onEnter"), this);

    }

    mLastExecutedIndex = mCurrentIndex;

    mCurrentIndex = indexmin;

    if (mCurrentIndex>-1) {

        int to_sleep = tmin - absoluteAge();

        if (to_sleep>0)

            sleep(to_sleep);

    }

    mScheduledHarvest = 0.; // reset

    mLastExecution = ForestManagementEngine::instance()->currentYear();

    return mCurrentIndex > -1;

}

void FMStand::notifyTreeRemoval(Tree *tree, int reason)

{

    double removed_volume = tree->volume();

    mVolume -= removed_volume/area();

    // for MAI calculations: store removal regardless of the reason

    mRemovedVolumeDecade+=removed_volume / area();

    mRemovedVolumeTotal+=removed_volume / area();

    Tree::TreeRemovalType r = Tree::TreeRemovalType (reason);

    if (r == Tree::TreeDeath)

        return; // do nothing atm

    if (r==Tree::TreeHarvest) {

        // regular harvest

        if (currentActivity()) {

            if (currentFlags().isFinalHarvest())

                mFinalHarvested += removed_volume;

            else

                mThinningHarvest += removed_volume;

        }

    }

    if (r==Tree::TreeDisturbance) {

        // if we have an active salvage activity, then store

        mDisturbed += removed_volume;

        // check if we have an (active) salvage activity; both the activity flags and the stand flags need to be "enabled"

        if (mSTP->salvageActivity() && mSTP->salvageActivity()->standFlags().enabled() && mSTP->salvageActivity()->standFlags(this).enabled() ) {

            if (mSTP->salvageActivity()->evaluateRemove(tree)) {

                mFinalHarvested += removed_volume;

                tree->setIsHarvested(); // set the flag that the tree is removed from the forest

            }

        }

    }

}

bool FMStand::notifyBarkBeetleAttack(double generations, int infested_px_per_ha)

{

    // check if we have an (active) salvage activity; both the activity flags and the stand flags need to be "enabled"

    if (mSTP->salvageActivity() && mSTP->salvageActivity()->standFlags().enabled() && mSTP->salvageActivity()->standFlags(this).enabled()) {

        return mSTP->salvageActivity()->barkbeetleAttack(this, generations, infested_px_per_ha);

    }

    return false;

}

void FMStand::sleep(int years_to_sleep)

{

    mYearsToWait = qMax(mYearsToWait, qMax(years_to_sleep,0));

}

double FMStand::calculateMAI()

{

    // MAI: delta standing volume + removed volume, per year

    // removed volume: mortality, management, disturbances

    mMAIdecade = ((mVolume - mLastMAIVolume) + mRemovedVolumeDecade) / 10.;

    if (absoluteAge()>0)

        mMAItotal = (mVolume + mRemovedVolumeTotal) / absoluteAge();

    mLastMAIVolume = mVolume;

    // reset counters

    mRemovedVolumeDecade = 0.;

    return meanAnnualIncrementTotal();

}

double FMStand::basalArea(const QString &species_id) const

{

    foreach (const SSpeciesStand &sd, mSpeciesData)

        if (sd.species->id()==species_id)

            return sd.basalArea;

    return 0.;

}

double FMStand::relBasalArea(const QString &species_id) const

{

    foreach (const SSpeciesStand &sd, mSpeciesData)

        if (sd.species->id()==species_id)

            return sd.relBasalArea;

    return 0.;

}

void FMStand::setAbsoluteAge(const double age)

{

    mRotationStartYear = ForestManagementEngine::instance()->currentYear() - age;

    mAge = age;

}

// storage for properties (static)

QHash<const FMStand*, QHash<QString, QJSValue> > FMStand::mStandPropertyStorage;

void FMStand::setProperty(const QString &name, QJSValue value)

{

    // save a property value for the current stand

    mStandPropertyStorage[this][name] = value;

}

QJSValue FMStand::property(const QString &name) const

{

    // check if values are already stored for the current stand

    if (!mStandPropertyStorage.contains(this))

        return QJSValue();

    // check if something is stored for the property name (return a undefined value if not)

    if (!mStandPropertyStorage[this].contains(name))

        return QJSValue();

    return mStandPropertyStorage[this][name];

}

QStringList FMStand::info()

{

    QStringList lines;

    lines << QString("id: %1").arg(id())

          << QString("unit: %1").arg(unit()->id());

    lines  << "-" << unit()->info() << "/-"; // sub sections

    if (currentActivity()) {

        lines << QString("activity: %1").arg(currentActivity()->name()) << "-" << currentActivity()->info();

        // activity properties

        lines << QString("active: %1").arg(currentFlags().active());

        lines << QString("enabled: %1").arg(currentFlags().enabled());

        lines << QString("simulate: %1").arg(currentFlags().isDoSimulate());

        lines << QString("execute immediate: %1").arg(currentFlags().isExecuteImmediate());

        lines << QString("final harvest: %1").arg(currentFlags().isFinalHarvest());

        lines << QString("use scheduler: %1").arg(currentFlags().isScheduled());

        lines << QString("in scheduler: %1").arg(currentFlags().isPending());

        lines <<  "/-";

    }

    lines << QString("agent: %1").arg(unit()->agent()->type()->name());

    lines << QString("STP: %1").arg(stp()?stp()->name():QStringLiteral("-"));

    lines << QString("U (yrs): %1").arg(U());

    lines << QString("thinning int.: %1").arg(thinningIntensity());

    lines << QString("last update: %1").arg(lastUpdate());

    lines << QString("sleep (years): %1").arg(sleepYears());

    lines << QString("scheduled harvest: %1").arg(scheduledHarvest());

    lines << QString("basal area: %1").arg(basalArea());

    lines << QString("volume: %1").arg(volume());

    lines << QString("age: %1").arg(age());

    lines << QString("absolute age: %1").arg(absoluteAge());

    lines << QString("N/ha: %1").arg(stems());

    lines << QString("MAI (decadal) m3/ha*yr: %1").arg(meanAnnualIncrement());

    lines << "Basal area per species";

    for (int i=0;i<nspecies();++i) {

        lines << QString("%1: %2").arg(speciesData(i).species->id()).arg(speciesData(i).basalArea);

    }

    lines  << "All activities" << "-";

    for (QVector<ActivityFlags>::const_iterator it = mStandFlags.constBegin(); it!=mStandFlags.constEnd(); ++it) {

        const ActivityFlags &a = *it;

        lines << QString("%1 (active): %2").arg(a.activity()->name()).arg(a.active())

                 << QString("%1 (enabled): %2").arg(a.activity()->name()).arg(a.enabled());

    }

    lines << "/-";

    // stand properties

    if (mStandPropertyStorage.contains(this)) {

        QHash<QString, QJSValue> &props = mStandPropertyStorage[this];

        lines << QString("properties: %1").arg(props.size()) << "-";

        QHash<QString, QJSValue>::const_iterator i = props.constBegin();

        while (i != props.constEnd()) {

            lines << QString("%1: %2").arg(i.key()).arg(i.value().toString());

            ++i;

        }

        lines << "/-";

    }

    // scheduler info

    lines << unit()->constScheduler()->info(id());

    return lines;

}

void FMStand::newRotatation()

{

    mLastRotationAge = absoluteAge();

    mRotationStartYear = ForestManagementEngine::instance()->currentYear(); // reset stand age to 0.

    mRemovedVolumeTotal = 0.;

    mRemovedVolumeDecade = 0.;

    mLastMAIVolume = 0.;

    mMAIdecade = 0.;

    mMAItotal = 0.;

    // use default values

    setU( unit()->U() );

    setThinningIntensity( unit()->thinningIntensity() );

    unit()->agent()->type()->agentUpdateForStand(this, QString(), 0); // update at age 0? maybe switch to new STP?

}

SSpeciesStand &FMStand::speciesData(const Species *species)

{

    for (int i=0;i<mSpeciesData.count(); ++i)

        if (mSpeciesData[i].species == species)

            return mSpeciesData[i];

    mSpeciesData.append(SSpeciesStand());

    mSpeciesData.last().species = species;

    return mSpeciesData.last();

}

} // namespace