Rabbit_Population_Manager.cpp

Go to the documentation of this file.

/*
*******************************************************************************************************
 
Copyright (c) 2015, Christopher John Topping, Faarupvej 54, DK-8410 Rønde
ADAMA Makhteshim Ltd., PO Box 60, Beer-Sheva 84100, Israel
 
All rights reserved.
 
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
 
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of the FreeBSD Project.
 
********************************************************************************************************
 
*/
//---------------------------------------------------------------------------
 
//#define  __RABBITBREEDINGSEASONCHECK
//#define __RABBITDEBUG2
 
 
#include <iostream>
#include <fstream>
#include <vector>
#include <math.h>
#include "../BatchALMaSS/ALMaSS_Setup.h"
#include "../ALMaSSDefines.h"
#include "../Landscape/ls.h"
#include "../BatchALMaSS/PopulationManager.h"
#include "../BatchALMaSS/AOR_Probe.h"
#include "../Rabbit/Rabbit.h"
#include "../Rabbit/Rabbit_Population_Manager.h"
#include "../BatchALMaSS/BoostRandomGenerators.h"
//---------------------------------------------------------------------------------------
 
extern boost::variate_generator<base_generator_type&, boost::uniform_real<> > g_rand_uni;
extern CfgFloat cfg_maxForageHeight;
extern CfgFloat cfg_dispersalmortperm;
extern CfgFloat cfg_rabbitminkits;
extern CfgFloat cfg_rabbitmaxkits;
extern CfgFloat cfg_minForageDigestability;
extern CfgBool cfg_RipleysOutput_used;
extern CfgBool cfg_AOROutput_used;
extern CfgFloat cfg_rabbit_pesticidedegradationrate;
//---------------------------------------------------------------------------------------
 
static CfgInt cfg_RabbitStartNos("RABBIT_STARTNOS",CFG_CUSTOM,25000);
CfgInt cfg_warrenfixedsize("RABBIT_WARRENFIXEDSIZE",CFG_CUSTOM,45);
static CfgBool cfg_warrenreadlocations("RABBIT_WARRENREADLOCATIONS", CFG_CUSTOM, true);
static CfgStr cfg_warrenlocationsfile("RABBIT_WARRENLOCATIONSFILE", CFG_CUSTOM, "RabbitWarrenLocations.txt");
static CfgInt cfg_warrenoccupancyrecordday("RABBIT_WARRENOCCUPANCYRECORDDAY", CFG_CUSTOM, 60);
static CfgInt cfg_maxWarrenNetworkDist("RABBIT_MAXWARRENNETWORKDIST",CFG_CUSTOM,1500);
static CfgFloat cfg_rabbitminbreedingtemp("RABBIT_MINBREEDINGTEMP",CFG_CUSTOM,300.0);
static CfgBool cfg_rabbit_use_fixed_soiltype("RABBIT_USE_FIXED_SOIL_TYPE", CFG_CUSTOM, true); 
static CfgInt cfg_rabbit_fix_soiltype("RABBIT_FIX_SOIL_TYPE", CFG_CUSTOM, 0); // 0 chalk, 1 sand
CfgInt cfg_rabbitdensitydependencedelay( "RABBIT_DENDEPPERIOD", CFG_CUSTOM, 53 );
CfgInt cfg_rabbitdiseasedensitydependencedelay( "RABBITDISEASEDENDEPPERIOD", CFG_CUSTOM, 90 );
static CfgFloat cfg_rabbitdailygrowthparam1("RABBIT_DAILYGROWTHPARAM_ONE", CFG_CUSTOM, 1127.616084);
static CfgFloat cfg_rabbitdailygrowthparam2("RABBIT_DAILYGROWTHPARAM_TWO", CFG_CUSTOM, -0.013143202);
static CfgFloat cfg_rabbitdailygrowthparam3("RABBIT_DAILYGROWTHPARAM_THREE", CFG_CUSTOM, 0.0);
static CfgFloat cfg_rabbitminimumforagetemp( "RABBIT_MINFORAGETEMP", CFG_CUSTOM, -0.83 );
static CfgFloat cfg_rabbitmaxforagerainfall( "RABBIT_MAXFORAGERAINFALL", CFG_CUSTOM, 3.6 );
CfgBool cfg_RabbitUseNatalDispersalRecord( "RABBIT_USENATALDISPERSALRECORD", CFG_CUSTOM, false );
CfgInt cfg_RabbitUseNatalDispersalRecordAge( "RABBIT_USENATALDISPERSALRECORDAGE", CFG_CUSTOM, 30*6 );
CfgBool cfg_RabbitUseReproOutput("RABBIT_USEREPROOUTPUT", CFG_CUSTOM, false);
// Globals created from input parameters - this is just to avoid unnecessary calculations later. NB these must be initialized in the code after program start 
int g_land_width = 0;
 
//***************************************************************************
//*************** RABBIT POPULATION MANAGER CODE ****************************
//***************************************************************************
 
Rabbit_Population_Manager::Rabbit_Population_Manager(Landscape* L) : Population_Manager(L)
{
    m_ListNames[0] = "Young";
    m_ListNames[1] = "Juvenile";
    m_ListNames[2] = "Male";
    m_ListNames[3] = "Female";
    m_ListNames[4] = "Possible Warrens";
    m_ListNameLength = 5;
    g_land_width = m_TheLandscape->SupplySimAreaWidth(); // Needs the landscape to be square!!!
 
    strcpy( m_SimulationName, "Rabbit Simulation" );
    struct_Rabbit* sp = new struct_Rabbit;
    sp->m_NPM = this;
    sp->m_L = m_TheLandscape;
    sp->m_age = 100+random(265);
    sp->m_Warren = NULL;
    sp->m_x2 = -1;
    sp->m_y2 = -1;
    sp->m_weightage = 300; // Should not be weak animals to start with
 
    for (int i = 0; i< cfg_RabbitStartNos.value(); i++)
    {
        sp->m_x = random(SimW);
        sp->m_y = random(SimH);
        CreateObjects(rob_Male,NULL,sp,1); // 
    }
    for (int i=0; i< cfg_RabbitStartNos.value(); i++) 
    {
        sp->m_x = random(SimW);
        sp->m_y = random(SimH);
        CreateObjects(rob_Female,NULL,sp,1); // 
    }
    delete sp;
    // Load parameters
    Rabbit_Warren::m_maxForageHeight = cfg_maxForageHeight.value();
    Rabbit_Base::m_dispersalmortperm = cfg_dispersalmortperm.value();
    m_warrenfixedsize = cfg_warrenfixedsize.value();
    m_warrenfixedsizediv2 = m_warrenfixedsize / 2;
    m_reproswitchbuffer = 30;
    m_rabbitBreedingSeason = false;
 
    if (cfg_warrenreadlocations.value())
    {
        LoadWarrenLocations();
    }
    else
    {   
        PreProcessWarrenLocations();
        SaveWarrenLocations();
    }
    CreateLocalWarrenNetworkLists();
    // Open output files
    if (cfg_RipleysOutput_used.value()) {
        OpenTheRipleysOutputProbe("");
    }
    PesticideDeathRecordOutputOpen();
    if (cfg_RabbitUseReproOutput.value()) ReproOutputRecordOutputOpen();
    if (cfg_RabbitUseNatalDispersalRecord.value()) NatalDispersalRecordOutputOpen();
    WarrenOccupancyRecordOutputOpen();
    // Assign static variable
    AssignStaticVariables();
    // Create the rabbit growth rate lookup data
    for (int i = 1; i < 300; i++) {
        m_RabbitGrowth[i] = cfg_rabbitdailygrowthparam1.value() - (cfg_rabbitdailygrowthparam1.value()*exp(cfg_rabbitdailygrowthparam2.value()*i)) + cfg_rabbitdailygrowthparam3.value();
    }
    // Now fill this in for up to 10 years (no rabbit should ever live that long )
    for (int i = 300; i < 3350; i++) {
        m_RabbitGrowth[ i ] = cfg_rabbitdailygrowthparam1.value() - (cfg_rabbitdailygrowthparam1.value()*exp( cfg_rabbitdailygrowthparam2.value()*300 )) + cfg_rabbitdailygrowthparam3.value();
    }
 
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::AssignStaticVariables(void)
{
    struct_Rabbit* sp = new struct_Rabbit;
    sp->m_NPM = this;
    sp->m_L = m_TheLandscape;
    sp->m_age = 100 + random(265);
    sp->m_Warren = NULL;
    Rabbit_Female* rf = new Rabbit_Female(0,0,-1,-1,m_TheLandscape,this,0, 300,NULL);
    Rabbit_Warren* rw = new Rabbit_Warren(0, 0, m_TheLandscape, this, 0, NULL);
    rw->m_maxForageHeight = cfg_maxForageHeight.value();
    rw->m_minForageDigestability = cfg_minForageDigestability.value();
    rf->m_dispersalmortperm = cfg_dispersalmortperm.value();
    rf->SetMinKits(cfg_rabbitminkits.value());
    rf->SetMaxKits(cfg_rabbitmaxkits.value() - cfg_rabbitminkits.value()); // Here we actually use the difference to calculate range
    rf->m_pesticidedegradationrate = cfg_rabbit_pesticidedegradationrate.value(); // default of 0.0 will remove all body burden pesticide at the end of each day
    delete rw;
    delete rf;
}
//---------------------------------------------------------------------------
 
Rabbit_Population_Manager::~Rabbit_Population_Manager(void)
{
    PesticideDeathRecordOutputClose();
    WarrenOccupancyRecordOutputClose();
    if (cfg_RabbitUseReproOutput.value()) {
        ReproOutputRecordOutputClose();
        LifetimeReproAnalysis();
    }
    if (cfg_RabbitUseNatalDispersalRecord.value()) {
        NatalDispersalRecordOutputClose();
        NatalDispersalAnalysis();
    }
    WarrenOutputAnalysis();
    for (unsigned i = 0; i < TheArray.size(); i++) {
    for ( unsigned j = 0; j < TheArray[ i ].size(); j++ ) {
        // This just stops a DEBUG error
        TheArray[ i ] [ j ]->KillThis();
    }
  }
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::CreateObjects(RabbitObjectTypes ob_type, TAnimal * pvo, struct_Rabbit * a_data, int a_number)
{
    switch (ob_type)
    {
    case rob_Young:
       Rabbit_Young*  new_RabbitY;
       for (int i=0; i<a_number; i++)
       {
           new_RabbitY = new Rabbit_Young(a_data->m_x, a_data->m_y, a_data->m_x2, a_data->m_y2, dynamic_cast<Rabbit_Female*>(pvo), a_data->m_L, a_data->m_NPM, a_data->m_Warren);
           TheArray[ob_type].push_back(new_RabbitY);
           dynamic_cast<Rabbit_Female*>(pvo)->AddYoung(new_RabbitY);
       }
       break;
    case rob_Juvenile:
       Rabbit_Juvenile*  new_RabbitJ;
       for (int i=0; i<a_number; i++)
       {
           new_RabbitJ = new Rabbit_Juvenile(a_data->m_x, a_data->m_y, a_data->m_x2, a_data->m_y2, dynamic_cast<Rabbit_Female*>(pvo), a_data->m_L, a_data->m_NPM, a_data->m_age, a_data->m_weightage, a_data->m_Warren);
           TheArray[ob_type].push_back(new_RabbitJ);
       }
       break;
    case rob_Male:
       Rabbit_Male*  new_RabbitM;
       for (int i=0; i<a_number; i++)
       {
           new_RabbitM = new Rabbit_Male(a_data->m_x, a_data->m_y, a_data->m_x2, a_data->m_y2, a_data->m_L, a_data->m_NPM, a_data->m_age, a_data->m_weightage, a_data->m_Warren);
           TheArray[ob_type].push_back(new_RabbitM);
       }
       break;
    case rob_Female:
       Rabbit_Female*  new_RabbitF;
       for (int i=0; i<a_number; i++)
       {
           new_RabbitF = new Rabbit_Female(a_data->m_x, a_data->m_y, a_data->m_x2, a_data->m_y2, a_data->m_L, a_data->m_NPM, a_data->m_age, a_data->m_weightage, a_data->m_Warren);
           TheArray[ob_type].push_back(new_RabbitF);
       }
       break;
    case rob_Warren:
       Rabbit_Warren*  new_RabbitW;
       for (int i=0; i<a_number; i++)
       {
           new_RabbitW = new Rabbit_Warren(a_data->m_x, a_data->m_y, a_data->m_L, a_data->m_NPM, m_warrenfixedsize, a_data->m_soil); 
           if (new_RabbitW->GetCarryingCapacity()>0)   TheArray[ob_type].push_back(new_RabbitW);
           else delete new_RabbitW;
       }
       break;
    }
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::DoFirst()
{
#ifdef __RABBITDEBUG2
    if (m_TheLandscape->SupplyYearNumber()== 1000)
        if (m_TheLandscape->SupplyDayInYear() == 1) {
            struct_Rabbit* sp = new struct_Rabbit;
            sp->m_NPM = this;
            sp->m_L = m_TheLandscape;
            sp->m_age = 100 + random( 265 );
            sp->m_Warren = NULL;
            sp->m_x2 = -1;
            sp->m_y2 = -1;
            sp->m_weightage = 300; // Should not be weak animals to start with
 
            for (int i = 0; i< 500; i++) {
                sp->m_x = random( SimW );
                sp->m_y = random( SimH );
                CreateObjects( rob_Male, NULL, sp, 1 ); // 
            }
            for (int i = 0; i< 500; i++) {
                sp->m_x = random( SimW );
                sp->m_y = random( SimH );
                CreateObjects( rob_Female, NULL, sp, 1 ); // 
            }
            delete sp;
        }
#endif
    CheckForRabbitBreedingConditions(); 
    for (int r = (int)rob_Young; r < (int)rob_foobar; r++)
    {
        m_PesticideDeaths[r] = 0;
    }
    // Manage warren disease
    if (m_TheLandscape->SupplyDayInYear() % cfg_rabbitdiseasedensitydependencedelay.value() == 0)
    {
        for (unsigned i = 0; i < TheArray[rob_Warren].size(); i++) dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i])->CalcDisease();
    }
    // Determine whether today is forage day
    double rainfall = m_TheLandscape->SupplyRain();
    double temp = m_TheLandscape->SupplyTemp();
    if ((temp<cfg_rabbitminimumforagetemp.value()) || (rainfall>cfg_rabbitmaxforagerainfall.value())) m_forageday = false; else m_forageday = true;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::DoLast()
{
    int today = m_TheLandscape->SupplyDayInYear();
    // Do the output jobs
    // Pesticide deaths
    PesticideDeathRecordOutput(); 
    // Warren occupancy
    if (today == cfg_warrenoccupancyrecordday.value()) WarrenOccupancyRecordOutput();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::PreProcessWarrenLocations()
{
    cout << "Preprocessing Warren Locations" << endl;
    int x = 0;
    int y = 0;
    int soil = 0;
    for (int j=y; j<SimH-m_warrenfixedsize; j++)
    {
        for (int i=x; i<SimW-m_warrenfixedsize; i++)
        {
            if ( WarrenLegalPos(i,j) )
            {
                if (cfg_rabbit_use_fixed_soiltype.value()) soil = cfg_rabbit_fix_soiltype.value();
                else
                {
                    soil = m_TheLandscape->SupplySoilTypeR(i, j);
                }
                double pct = AssessPctForage(i,j);
                if ((soil) != 3)  // Soil type 3 is unsuitable therefore skipped
                {
                    // Soil type 1 is sandy, 0 is chalk or otherwise hard soil
                    if (pct >= 1.0)
                    {
                        struct_Rabbit sR;
                        sR.m_NPM = this;
                        sR.m_L = m_TheLandscape;
                        sR.m_x = i;
                        sR.m_y = j;
                        sR.m_soil = soil;
                        CreateObjects(rob_Warren, NULL, &sR, 1);
                        --i += m_warrenfixedsize;
                        if (static_cast<unsigned>(TheArray.size()) % 1000 == 0) cout << TheArray.size() << '\t';
                    }
                }
            }
        }
    }
}
//---------------------------------------------------------------------------
 
bool Rabbit_Population_Manager::WarrenLegalPos(int &a_x, int a_y)
{
    if (ClassifyHabitat(m_TheLandscape->SupplyElementType(a_x+m_warrenfixedsizediv2,a_y+m_warrenfixedsizediv2)) != torh_Forage) return false;
    // If it is not in a warren already then we are legal here. 
    int sz = (int) TheArray[rob_Warren].size();
    for (int i=0; i< sz; i++) 
    {
        APoint pt = dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i])->SupplyPoint();   
        if ((abs(pt.m_x-a_x) < m_warrenfixedsize) && (abs(pt.m_y-a_y) < m_warrenfixedsize))
        {
            a_x += m_warrenfixedsize;
            return false;
        }
    }
    return true;
}
//---------------------------------------------------------------------------
 
TTypesOfRabbitHabitat Rabbit_Population_Manager::ClassifyHabitat(TTypesOfLandscapeElement a_tole)
{
    static char error_num[ 20 ];
    switch (a_tole)
    {
        case tole_UnsprayedFieldMargin://31
        case tole_Field://20&30
        case tole_RoadsideVerge://13
        case tole_BeetleBank://141
            return torh_TemporaryForage;
            break;
        case tole_Railway://118
        case tole_FieldBoundary://160
        case tole_PermPastureLowYield://26
        case tole_PermPastureTussocky://27
        case tole_PermanentSetaside://33
        case tole_PermPasture://35
        case tole_NaturalGrassDry://110//case tole11
        case tole_PitDisused://75
        case tole_YoungForest://55
        case tole_HedgeBank://140
        case tole_Heath://94
        case tole_Parkland://14
        case tole_SandDune://101
        case tole_Orchard://56
        case tole_RoadsideSlope:  //201
        case tole_HeritageSite:  //208
        case tole_UnknownGrass:
        case  tole_Vildtager:
        case tole_MownGrass:
        case tole_Wasteland:
            return torh_Forage;
            break;
        case tole_Hedges://130
        case tole_Copse://41
        case tole_Scrub://70
        case tole_RiversideTrees://97
        case tole_DeciduousForest://40
        case tole_MixedForest://60
        case tole_ConiferousForest://50
        case tole_WoodlandMargin:
        case tole_WoodyEnergyCrop:
            return torh_Cover;
            break;
        case tole_UrbanPark:  //17
        case tole_NaturalGrassWet:
        case tole_Churchyard: //204
        case tole_Saltmarsh:  //206
        case tole_Stream:  //207
        case tole_Carpark:  //203
        case tole_IndividualTree:  //42
        case tole_MetalledPath:  //202
        case tole_StoneWall://15
        case tole_Fence: //225
        case tole_Garden://11//case tole20
        case tole_Track://123
        case tole_SmallRoad://122
        case tole_LargeRoad://121
        case tole_Building://5
        case tole_ActivePit://115
        case tole_Pond:
        case tole_Freshwater://90
        case tole_FishFarm:
        case tole_River://96
        case tole_Saltwater://80
        case tole_Coast://100
        case tole_UrbanNoVeg://8
        case tole_BuiltUpWithParkland://16
        case tole_AmenityGrass://12
        case tole_OrchardBand://57
        case tole_BareRock://59
        case tole_Marsh://95
        case tole_RiversidePlants://98
        case tole_PlantNursery:
        case tole_Pylon:
        case tole_WindTurbine:
        case  tole_PermPastureTussockyWet:
            return torh_Other;
            break;
        default:
            sprintf( error_num, "%d", a_tole );
            m_TheLandscape->Warn("Rabbit_Warren::InitEvaluation Unknown tole type: ", error_num);
            exit(1);
            break;
    }
}
//---------------------------------------------------------------------------
 
double Rabbit_Population_Manager::AssessPctForage(int a_x, int a_y)
{   
    int total = 0;
    int add = 0;
    int oldpoly = -1;
    for (int i=a_x; i<a_x + m_warrenfixedsize; i++)
    {
        for (int j=a_y; j<a_y + m_warrenfixedsize; j++)
        {
            int poly = m_TheLandscape->SupplyPolyRef(i, j);
            if (poly != oldpoly)
            {
                oldpoly = poly;
                TTypesOfLandscapeElement tole = m_TheLandscape->SupplyElementType(i, j);
                if (ClassifyHabitat(tole) == torh_Forage) add = 1; else add = 0;
            }
            total += add;
        }
    }
    double warrenarea = m_warrenfixedsize * m_warrenfixedsize;
    return (double) total/warrenarea;
}
 
void Rabbit_Population_Manager::SaveWarrenLocations()
{
    ofstream ofile(cfg_warrenlocationsfile.value());
    if ( !ofile.is_open() ) 
    {
        g_msg->Warn("Rabbit_Population_Manager::SaveWarrenLocations() Cannot open file: ", cfg_warrenlocationsfile.value());
      exit( 1 );
    }
    int sz = (int) TheArray[rob_Warren].size();
    ofile << sz << endl;
    for (int i=0; i< sz; i++) 
    {
        Rabbit_Warren* awarren = dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i]);
        ofile << awarren->Supply_m_Location_x() - m_warrenfixedsizediv2 << '\t' << awarren->Supply_m_Location_y() - m_warrenfixedsizediv2 << '\t' << awarren->GetSoilType() << endl;
    }
    ofile.close();
}
 
void Rabbit_Population_Manager::LoadWarrenLocations()
{
    ifstream ifile(cfg_warrenlocationsfile.value());
    if ( !ifile.is_open() ) 
    {
        g_msg->Warn("Rabbit_Population_Manager::LoadWarrenLocations() Cannot open file: ", cfg_warrenlocationsfile.value());
      exit( 1 );
    }
    int sz;
    ifile >> sz;
    for (int i=0; i< sz; i++) 
    {
        int x,y,s;
        ifile >> x >> y >> s;
        struct_Rabbit sR;
        sR.m_NPM = this;
        sR.m_L = m_TheLandscape;
        sR.m_x = x;
        sR.m_y = y;
        sR.m_soil = s;
        CreateObjects(rob_Warren,NULL,&sR,1); 
    }
    ifile.close();
}
 
 
void Rabbit_Population_Manager::CreateLocalWarrenNetworkLists( void )
{
    int sz = (int) TheArray[rob_Warren].size();
    for (int i=0; i< sz; i++) 
    {
        int closest = -1;
        int dist = SimH;
        int mDist = SimH;
        int ourX = TheArray[rob_Warren][i]->Supply_m_Location_x();
        int ourY = TheArray[rob_Warren][i]->Supply_m_Location_y();
        for (int j=0; j< sz; j++) 
        {
            if (i!=j) // Don't add ourselves
            {
                int theirX = TheArray[rob_Warren][j]->Supply_m_Location_x();
                int theirY = TheArray[rob_Warren][j]->Supply_m_Location_y();
                dist = (int) sqrt( (double)(((theirX-ourX)*(theirX-ourX)) + ((theirY-ourY)*(theirY - ourY))) );
                if (dist < cfg_maxWarrenNetworkDist.value()) 
                {
                    // Add to the network for this warren
                    LocalWarrenNewtorkEntry LWNE;
                    LWNE.m_aWarren = dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][j]);
                    LWNE.m_dist = dist;
                    dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i])->AddNetworkConnection( LWNE );
                }
                if (dist < mDist) 
                { 
                    mDist = dist; 
                    closest = j; 
                }
            }
        }
        if (dist > cfg_maxWarrenNetworkDist.value())
        {
            // Add the closest warren since there are no others
            LocalWarrenNewtorkEntry LWNE;
            LWNE.m_aWarren = dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][closest]);
            LWNE.m_dist = dist;
            dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i])->AddNetworkConnection( LWNE );
        }
    }
    for (int i=0; i< sz; i++) 
    {
        dynamic_cast<Rabbit_Warren*>(TheArray[rob_Warren][i])->NetworkEvaluation();
    }
}
 
 
Rabbit_Warren* Rabbit_Population_Manager::FindClosestWarren(int a_x, int a_y, int a_rank)
{
    int results[10];
    int dist;
    TAnimal* ptrs[10];
    for (int r = 0; r < 10; r++) {
        results[ r ] = SimH*SimW;
        ptrs[ r ] = NULL;
    }
    int sz = (int) TheArray[rob_Warren].size();
    for (int i=0; i< sz; i++) 
    {
        int X = TheArray[rob_Warren][i]->Supply_m_Location_x();
        int Y = TheArray[rob_Warren][i]->Supply_m_Location_y();
        dist = (int)sqrt( ((double( X - a_x )*double( X - a_x )) + (double( Y - a_y )*double( Y - a_y ))) );
        for (int r=0; r<10; r++)
        {
            if (dist<=results[r])
            {
                for (int t=r+1; t<10; t++)
                {
                    results[t]=results[t-1];
                    ptrs[t]=ptrs[t-1];
                }
                results[r]=dist;
                ptrs[r]=TheArray[rob_Warren][i];
                break;
            }
        }
    }
    Rabbit_Warren*  rw = dynamic_cast<Rabbit_Warren*>(ptrs[ a_rank ]);
    return rw;
}
 
void Rabbit_Population_Manager::CheckForRabbitBreedingConditions( void )
{
    /*
    * Once the temperature comes up over a certain trigger the breeding seasons starts. Rapid switching is prevented by not making the test again for 120 days.
    * When the temperature in the future drops down again then the breeding season stops
    */
    if (--m_reproswitchbuffer > 0) return;
 
    if (!m_rabbitBreedingSeason)
    {
        if (m_TheLandscape->SupplyTempPeriod(g_date->Date()+60,60) > cfg_rabbitminbreedingtemp.value())
        {
            m_reproswitchbuffer = 120;
            m_rabbitBreedingSeason = true;
#ifdef __RABBITBREEDINGSEASONCHECK
            cout << g_date->GetYearNumber() << " " << g_date->DayInYear() << " " << m_TheLandscape->SupplyTempPeriod( g_date->Date()+60, 60 ) << " - " << "Started" << endl;
#endif
        }
    }
    else
    {
if (m_TheLandscape->SupplyTempPeriod(g_date->Date()+60,60) < cfg_rabbitminbreedingtemp.value() ) 
        {
            m_reproswitchbuffer = 90;
            m_rabbitBreedingSeason = false;
#ifdef __RABBITBREEDINGSEASONCHECK
            cout << g_date->GetYearNumber() << " " << g_date->DayInYear() << " " << m_TheLandscape->SupplyTempPeriod( g_date->Date()+30, 90 ) << " - " << "Stopped" << endl;
#endif
        }
    }
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::PesticideDeathRecord(RabbitObjectTypes a_ob)
{
    m_PesticideDeaths[a_ob]++;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::WarrenOccupancyRecordOutput() {
    ofstream warrenrecord( "WarrenOccupancyRecord.txt", ios::app );
    double occupancy = 0;
    double var_occupancy = 0;
    double sumsqr = 0.0;
    double totalrabbits = 0;
    double burrows = 0;
    double disease = 0;
    int warrensoccupied = 0;
    int kits = 0;
    int juvs = 0;
    int ads = 0;
    int breedingfemales = 0;
    double breedingfemalesAv = 0;
    int oneyroldfemales = 0;
    int nonbreedingfemales = 0;
    int litters = 0;
    double lrc = 0;
    double n = (double)TheArray[ rob_Warren ].size();
    for (unsigned i = 0; i < TheArray[ rob_Warren ].size(); i++) {
        Rabbit_Warren* warren = dynamic_cast<Rabbit_Warren*>(TheArray[ rob_Warren ][ i ]);
        double occ = warren->GetCarryingCapacityRatio();
        occupancy += occ;
        sumsqr += occ * occ;
        totalrabbits += warren->GetPopulationSize();
        if (warren->GetPopulationSize()>0) warrensoccupied++;
        burrows += warren->GetCarryingCapacity();
        disease += warren->GetDiseaseConstant();
        kits += warren->GetRabbitProductionRecord( rob_Young );
        juvs += warren->GetRabbitProductionRecord( rob_Juvenile );
        ads += warren->GetRabbitProductionRecord( rob_Male );
        ads += warren->GetRabbitProductionRecord( rob_Female );
        warren->UpdateThisYearsBreeders();
        breedingfemales += warren->GetThisYearsBreeders();
        breedingfemalesAv += warren->GetThisYearsBreedersAv();
        oneyroldfemales += warren->GetThisYears1yrOldFemales();
        nonbreedingfemales += warren->GetThisYearsNonBreeders();
        litters += warren->GetLittersThisYear();
        lrc += warren->GetLitterReabsortionConst();
        warren->ResetAllRabbitProductionRecord();
    }
    double littersperfemale = 0;
    double CV_occ = 0;
    if ((occupancy > 0) && (ads > 0)) {
        var_occupancy = (sumsqr - (occupancy*occupancy) / n) / (n - 1);
        if (var_occupancy < 0) var_occupancy = 0.0; // Rounding errors cause this
        occupancy /= n;
        burrows /= n;
        disease /= n;
        totalrabbits /= n;
        kits /= (int)n;
        juvs /= (int)n;
        ads /= (int)n;
        CV_occ = sqrt( var_occupancy ) / occupancy;
        if (breedingfemales > 0) littersperfemale = litters / double( breedingfemales );
    }
    else {
        var_occupancy = 0;
        occupancy = 0;
        burrows = 0;
        disease = 0;
        totalrabbits = 0;
        kits = 0;
        juvs = 0;
        ads = 0;
        CV_occ = 0;
    }
    lrc /= n;
    warrenrecord << m_TheLandscape->SupplyYearNumber() << '\t' << m_TheLandscape->SupplyDayInYear() << '\t' << (int)n << '\t' << warrensoccupied << '\t' << burrows << '\t' << totalrabbits << '\t' << disease << '\t' << occupancy << '\t' << CV_occ << '\t' << kits << '\t' << juvs << '\t' << ads;
    for (int r = (int)rob_Young; r < rob_Warren; r++) warrenrecord << '\t' << TheArray[ r ].size();
    warrenrecord << '\t' << breedingfemales << '\t' << breedingfemalesAv << '\t' << oneyroldfemales << '\t' << nonbreedingfemales << '\t' << littersperfemale << '\t' << lrc << '\t' << litters << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::PesticideDeathRecordOutput()
{
    m_PesticideDeathOFile << m_TheLandscape->SupplyYear() << '\t' << m_TheLandscape->SupplyDayInYear() << '\t';
    for (int r = (int)rob_Young; r < (int)rob_foobar; r++)
    {
        m_PesticideDeathOFile << '\t' << m_PesticideDeaths[r];
    }
    m_PesticideDeathOFile << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::PesticideDeathRecordOutputOpen()
{
    m_PesticideDeathOFile.open("RabbitPesticideDeathRecord.txt", ios::out);
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::WarrenOccupancyRecordOutputOpen()
{
    m_WarrenOccupancyFile.open("WarrenOccupancyRecord.txt", ios::out);
    m_WarrenOccupancyFile << "Year" << '\t' << "Day" << '\t' << "No. Warrens" << '\t' << "No.Occ.Warrens" << '\t' << "Mean Burrows" << '\t' << "Mean Rabbits" << '\t' << "Mean Disease" 
        << '\t' << "Mean Occupancy" << '\t' << "CV Occupancy" << '\t' << "Mean Kits" << '\t' << "Mean Juvs" << '\t' << "Mean Ads" << '\t' << "TotalYoung" << '\t' << "TotalJuvs" 
        << '\t' << "TotalMales" << '\t' << "TotalFemales" << '\t' << "BreedingFemales" << '\t' << "BreedingDensityAv" << '\t' << "1yr old females" << '\t' << "NonBreedingFems.>1yr" 
        << '\t' << "Littersperfemale" << '\t' << "lrc" << '\t' << "litters" << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::PesticideDeathRecordOutputClose()
{
    m_PesticideDeathOFile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::ReproOutputRecordOutput(Rabbit_Female* a_female)
{
    m_ReproOutputFile << m_TheLandscape->SupplyYear() << '\t' << m_TheLandscape->SupplyDayInYear() << '\t';
    m_ReproOutputFile << '\t' << a_female->GetAge() << '\t' << a_female->GetTotalOffspring() << '\t' << a_female->GetTotalLitters() << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::ReproOutputRecordOutputOpen()
{
    m_ReproOutputFile.open( "RabbitLifetimeReproRecord.txt", ios::out );
    m_ReproOutputFile << "Year" << '\t' << "Day" << '\t';
    m_ReproOutputFile << '\t' << "Age" << '\t' << "Offpring" << '\t' << "Litters" << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::ReproOutputRecordOutputClose()
{
    m_ReproOutputFile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::NatalDispersalRecordOutput(Rabbit_Base* a_rabbit)
{
    m_NatalDispersalFile << m_TheLandscape->SupplyYear() << '\t' << m_TheLandscape->SupplyDayInYear() << '\t';
    APoint pt = a_rabbit->GetBornLocation();
    m_NatalDispersalFile << pt.m_x << '\t' << pt.m_y << '\t' << a_rabbit->Supply_m_Location_x() << '\t' << a_rabbit->Supply_m_Location_y() << '\t' << a_rabbit->GetRabbitType() << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::NatalDispersalRecordOutputOpen()
{
    m_NatalDispersalFile.open("RabbitNatalDispersalRecord.txt", ios::out);
    m_NatalDispersalFile << "Year" << '\t' << "Day" << '\t';
    m_NatalDispersalFile << "Born_X" << '\t' << "Born_Y" << '\t' << "New_X" << '\t' << "New_Y" << '\t' << "Rabbit_Type" << endl;
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::NatalDispersalRecordOutputClose()
{
    m_NatalDispersalFile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::WarrenOccupancyRecordOutputClose()
{
    m_WarrenOccupancyFile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::WarrenOutputAnalysis() {
    double Year, Day, NoWarrens, NoOccWarrens, MeanBurrows, MeanRabbits, MeanDisease, MeanOccupancy, CVOccupancy, MeanKits, MeanJuvs, MeanAds, TotalYoung, TotalJuvs, TotalMales,
        TotalFemales, BreedingFemales, yroldfemales, NonBreedingFems;
    vector<double>  totalfemales;
    vector<double>  totalbrfemalesdensity;
    vector<double>  nonbreedingdensity;
    vector<double>  yrolddensity;
    vector<double>  breedingdensity;
    vector<double>  yearlingpct;
    vector<double>  kitsperfemale;
 
    double meantotalbrpop;
    double maxbreedingpop = -1;
    double minbreedingpop = 100000;
    double maxyroldspct = -1;
    double minyroldspct = 100000;
    double correlation;
    double correlation_b;
    double littersperfemale;
    double lrc;
    double BreedingDensityAv;
    double litters;
 
 
    ifstream WarrenOccupancyFile( "WarrenOccupancyRecord.txt", ios::in );
    // Read the header line
    char st[ 2000 ];
    WarrenOccupancyFile.getline( st, 2000 );
    while (!WarrenOccupancyFile.eof()) {
        WarrenOccupancyFile >> Year >> Day >> NoWarrens >> NoOccWarrens >> MeanBurrows >> MeanRabbits >> MeanDisease
            >> MeanOccupancy >> CVOccupancy >> MeanKits >> MeanJuvs >> MeanAds >> TotalYoung >> TotalJuvs
            >> TotalMales >> TotalFemales >> BreedingFemales >> BreedingDensityAv >> yroldfemales >> NonBreedingFems >> littersperfemale >> lrc >> litters;
        // Here we create the actual data that we need
        // We want:
        // - correlation between total breeding pop and 1yr olds
        // - mean of total breeding pop
        // - mean of 1yr olds
        // - range max-min of total breeding pop
        // - range max-min of 1yr olds
        if ((yroldfemales + BreedingFemales + NonBreedingFems) > 0) {
            yearlingpct.push_back( yroldfemales / (yroldfemales + BreedingFemales + NonBreedingFems) );
            totalbrfemalesdensity.push_back( (yroldfemales + BreedingFemales + NonBreedingFems) / 2.0 );
            totalfemales.push_back( (BreedingFemales + NonBreedingFems) / 2.0 );
            breedingdensity.push_back( BreedingFemales / 2.0 );
            nonbreedingdensity.push_back( NonBreedingFems / 2.0 );
            yrolddensity.push_back( yroldfemales / 2.0 );
            kitsperfemale.push_back( littersperfemale );
        }
        else {
            // All dead so save the results
            ofstream ofile( "RabbitPOMSummary.txt", ios::app );
            ofile << 0
                << '\t' << 0
                << '\t' << 0
                << '\t' << 0
                << '\t' << 0
                << '\t' << 0
                << '\t' << 0 << endl;
            ofile.close();
            return;
        }
    }
    WarrenOccupancyFile.close();
    // Now we have all the useful data we can create the stats for this analysis
    // 
    // Get the ranges skipping the first 10 years
    int sz = static_cast<int>(yearlingpct.size());
    double sumx = 0;
    double sumx2 = 0;
    double sumy = 0;
    double sumy2 = 0;
    double sumxy = 0;
    double sumx_b = 0;
    double sumx2_b = 0;
    double sumy_b = 0;
    double sumy2_b = 0;
    double sumxy_b = 0;
    double n = sz - 1 - 10;
    // Do the correlation for yearlingpct and total female density
    for (int i = 10; i < sz - 1; i++) {
        if (totalbrfemalesdensity[ i ]>maxbreedingpop)maxbreedingpop = totalbrfemalesdensity[ i ];
        if (totalbrfemalesdensity[ i ] < minbreedingpop)minbreedingpop = totalbrfemalesdensity[ i ];
        if (yearlingpct[ i ] > maxyroldspct) maxyroldspct = yearlingpct[ i ];
        if (yearlingpct[ i ] < minyroldspct) minyroldspct = yearlingpct[ i ];
        sumx += totalbrfemalesdensity[ i ];
        sumx2 += totalbrfemalesdensity[ i ] * totalbrfemalesdensity[ i ];
        sumy += yearlingpct[ i ];
        sumy2 += yearlingpct[ i ] * yearlingpct[ i ];
        sumxy += totalbrfemalesdensity[ i ] * yearlingpct[ i ];
    }
    meantotalbrpop = sumx / n;
    double intermediate = (n*sumxy) - (sumx*sumy);
    correlation = intermediate / sqrt( ((n*sumx2) - (sumx*sumx))*((n*sumy2) - (sumy*sumy)) );
    //
    n = n - 1;
    for (int i = 10; i < sz - 2; i++) {
        // Here we need the numbers from the year before
        sumx_b += totalfemales[ i ];
        sumx2_b += totalfemales[ i ] * totalfemales[ i ];
        sumy_b += kitsperfemale[ i+1 ];
        sumy2_b += kitsperfemale[ i+1 ] * kitsperfemale[ i+1 ];
        sumxy_b += totalfemales[ i ] * kitsperfemale[ i+1 ];
    }
    intermediate = (n*sumxy_b) - (sumx_b*sumy_b);
    correlation_b = intermediate / sqrt( ((n*sumx2_b) - (sumx_b*sumx_b))*((n*sumy2_b) - (sumy_b*sumy_b)) );
    // Now save the results
    ofstream ofile( "RabbitPOMSummary.txt", ios::app );
    ofile << meantotalbrpop
          << '\t' << maxbreedingpop
          << '\t' << minbreedingpop
          << '\t' << maxyroldspct
          << '\t' << minyroldspct
          << '\t' << correlation
          << '\t' << correlation_b << endl;
    ofile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::LifetimeReproAnalysis() {
    /*
    % Success   0.66
    Average litters( of reproducers )   4.39
    Av.Total Offspring  21.13
    Mean Repro Age  348.26
    Max Repro Age   2069.00
    */
    int Year, Day, Age, Offspring, Litter;
    double NoSuccessfulReproducers = 0;
    double NoOffspring = 0;
    double TotalAdults = 0;
    double ZeroLitters = 0;
    double TotalLitters = 0;
    double AvNoYoungLitter = 0;
    double AvLittersPerFemale = 0;
    double AvTotalOffspring = 0;
    double pctSuccess = 0;
    double MeanReproAge = 0;
    double SumReproAge = 0;
    double MaxReproAge = 0;
    double sumunderage = 0;
    double sumunderageoffspring = 0;
 
    ifstream ReproFile( "RabbitLifetimeReproRecord.txt", ios::in );
    // Read the header line
    char st[ 2000 ];
    ReproFile.getline( st, 2000 );
    ReproFile >> Year >> Day >> Age >> Offspring >> Litter;
    while (!ReproFile.eof()) {
        if (Year > 2000) {
            if (Age >= 304) {
                TotalAdults++;
                SumReproAge += Age - 304;
                if (Litter > 0) {
                    TotalLitters += Litter;
                    NoSuccessfulReproducers++;
                    if (MaxReproAge < Age - 304) MaxReproAge = Age - 304;
                }
                else ZeroLitters++;
                NoOffspring += Offspring;
            }
            else if (Litter > 0) {
                sumunderage++;
                sumunderageoffspring += Offspring;
            }
        }
        ReproFile >> Year >> Day >> Age >> Offspring >> Litter;
    }
    if (TotalLitters > 0) AvNoYoungLitter = NoOffspring / TotalLitters;
    if (TotalAdults > 0) {
        MeanReproAge = SumReproAge / TotalAdults;
        pctSuccess = NoSuccessfulReproducers / TotalAdults;
        AvLittersPerFemale = TotalLitters / NoSuccessfulReproducers;
    }
    AvTotalOffspring = AvLittersPerFemale * AvNoYoungLitter;
    ofstream ofile( "RabbitPOMSummary.txt", ios::app );
    ofile << NoSuccessfulReproducers
        << '\t' << TotalLitters
        << '\t' << NoOffspring
        << '\t' << pctSuccess
        << '\t' << AvTotalOffspring
        << '\t' << AvNoYoungLitter
        << '\t' << AvLittersPerFemale
        << '\t' << MeanReproAge
        << '\t' << MaxReproAge
        << '\t' << sumunderage
        << '\t' << sumunderageoffspring << '\t';
    ofile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::NatalDispersalAnalysis() {
    int year, day, bornx, borny, newx, newy, sex;
    int sedmales = 0;
    int sedfemales = 0;
    int dispmales = 0;
    int dispfemales = 0;
    double malenataldisp = 0;
    double femalenataldisp = 0;
    ifstream DispFile( "RabbitNatalDispersalRecord.txt", ios::in );
    // Read the header line
    char st[ 2000 ];
    DispFile.getline( st, 2000 );
    DispFile >> year >> day >> bornx >> borny >> newx >> newy >> sex;
    while (!DispFile.eof()) {
        if (year > 2000) {
            if ((bornx != newx) || (borny != newy)) {
                if (sex == 2) dispmales++; else dispfemales++;
            }
            else if (sex == 2) sedmales++; else sedfemales++;
        }
        DispFile >> year >> day >> bornx >> borny >> newx >> newy >> sex;
    }
    if ((dispmales + sedmales)>0) malenataldisp = dispmales / double( dispmales + sedmales );
    if ((dispfemales + sedfemales)>0) femalenataldisp = dispfemales / double( dispfemales + sedfemales );
    ofstream ofile( "RabbitPOMSummary.txt", ios::app );
    ofile << malenataldisp  << '\t' << femalenataldisp << '\t';
    ofile.close();
}
//---------------------------------------------------------------------------
 
void Rabbit_Population_Manager::TheAOROutputProbe() {
    m_AOR_Probe->DoProbe(rob_Female);
}
//-----------------------------------------------------------------------------
 
// ***************************************************************************
// ************* END RABBIT POPULATION MANAGER CODE **************************
// ***************************************************************************