m_harvest.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file m_harvest.h
///
/// A nonlinear population model map with harvesting.
///
/// This map defines a nonlinear population model with harvesting.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// This file is part of my research software package. This 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 software 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 software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on February 9, 2008. Last revision: April 25, 2011.
 
 
#ifndef _CMGRAPHS_HARVESTING_H_
#define _CMGRAPHS_HARVESTING_H_
 
 
// some standard header files
#include <cmath>
 
// include local header files
#include "maptype.h"
#include "typeintv.h"
 
// for debugging
#include "chomp/system/textfile.h"
 
 
// --------------------------------------------------
// ------------- a population model map -------------
// --------------------------------------------------
 
/// This class defines a map for the nonlinear density dependent
/// overcompensatory Leslie population model.
/// The interpretation of parameters depends on the dimension.
class MapHarvesting: public MapType
{
public:
        /// The default constructor.
        MapHarvesting ();
 
        /// Computes the image of a box whose left and right coordinates
        /// are given. Fills in the images with the left and right
        /// coordinates of the image box.
        void compute (const double *xleft, const double *xright,
                double *yleft, double *yright, int dim,
                const spcCoord *coord, int subdiv) const;
 
}; /* class MapHarvesting */
 
// --------------------------------------------------
 
inline MapHarvesting::MapHarvesting ()
{
        testIntervals ();
        return;
} /* MapHarvesting::MapHarvesting */
 
inline void MapHarvesting::compute (const double *xleft, const double *xright,
        double *yleft, double *yright, int dim, const spcCoord *, int) const
{
        // this formula is programmed for phase space of dim 2 only
        if (dim != 2)
                throw "Population model with harvesting: dim = 2 only!";
 
        // determine the number of adults and juveniles
        IntervalType adults (xleft [0], xright [0]);
        IntervalType juveniles (xleft [1], xright [1]);
 
        // prepare the interval version of the number one
        IntervalType one (1.0, 1.0);
 
        // determine the parameters
        IntervalType harvestingRateA (getLeftParam (0), getRightParam (0));
        IntervalType harvestingRateJ (getLeftParam (1), getRightParam (1));
        IntervalType survivalRateA (getLeftParam (2), getRightParam (2));
        IntervalType survivalRateJ (getLeftParam (3), getRightParam (3));
        IntervalType alpha (getLeftParam (4), getRightParam (4));
        IntervalType beta (getLeftParam (5), getRightParam (5));
 
        // determine modifications of the model if any
        int model = static_cast<int> (getLeftParam (6));
        if (model & 0x001)
                harvestingRateJ = IntervalType (0.0, 0.0);
        else if (model & 0x002)
                harvestingRateA = IntervalType (0.0, 0.0);
        else if (model & 0x004)
                harvestingRateJ = harvestingRateA;
        if (model & 0x010)
                survivalRateJ = one - survivalRateA;
        else if (model & 0x020)
                survivalRateJ = survivalRateA;
        if (model & 0x100)
        {
                // use the cached value of beta if possible
                static IntervalType alphaSaved (0.0, 0.0);
                static IntervalType betaSaved (0.0, 0.0);
                if ((alphaSaved. leftBound () == alpha. leftBound ()) &&
                        (alphaSaved. rightBound () == alpha. rightBound ()))
                {
                        beta = betaSaved;
                }
 
                // compute the approximate value of beta if alpha is a point
                else if (alpha. leftBound () == alpha. rightBound ())
                {
                        static double e = M_E; //std::exp (1.0);
                        double b = alpha. leftBound () / e;
                        beta = IntervalType (b, b);
                        alphaSaved = alpha;
                        betaSaved = beta;
                }
 
                // compute the interval vector for beta otherwise
                else
                {
                        static IntervalType e = exp (one);
                        beta = alpha / e;
                        alphaSaved = alpha;
                        betaSaved = beta;
                }
        }
 
        IntervalType harvRemainingA (one - harvestingRateA);
        IntervalType harvRemainingJ (one - harvestingRateJ);
        IntervalType harvestedAdults (harvRemainingA * adults);
        IntervalType newJuveniles (alpha * harvestedAdults *
                exp (-beta * harvestedAdults));
        IntervalType newAdults
                (harvRemainingJ * survivalRateJ * juveniles +
                harvRemainingA * survivalRateA * adults);
 
        yleft [0] = newAdults. leftBound ();
        yright [0] = newAdults. rightBound ();
        yleft [1] = newJuveniles. leftBound ();
        yright [1] = newJuveniles. rightBound ();
 
//      chomp::homology::sbug << "Juveniles: " << juveniles <<
//              " -> " << newJuveniles << ".\n" << "Adults: " << adults <<
//              " -> " << newAdults << ".\n";
 
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapHarvesting::compute */
 
 
#endif // _CMGRAPHS_HARVESTING_H_