Go back to Richel Bilderbeek's homepage.

Go back to Richel Bilderbeek's C++ page.

 

 

 

 

 

(C++) QtPostProject

 

 

 

 

 

 

Technical facts

 

Application type(s)

Operating system(s) or programming environment(s)

IDE(s):

Project type:

C++ standard:

Compiler(s):

Libraries used:

 

 

 

 

 

Qt project file: ProjectChrisWileyDesktop.pro

 

#-------------------------------------------------
#
# Project created by QtCreator 2010-07-22T22:43:13
#
#-------------------------------------------------
SOURCES += \
    bird.cpp \
    matetally.cpp \
    speciestally.cpp \
    timepoint.cpp \
    timeseries.cpp \
    simulation.cpp \
    results.cpp \
    random.cpp \
    parameters.cpp \
    matingsystem.cpp \
    helperfunctions.cpp \
    femalesampling.cpp \
    densitydependentselection.cpp \
    main.cpp
HEADERS += \
    bird.h \
    matetally.h \
    speciestally.h \
    timepoint.h \
    timeseries.h \
    simulation.h \
    results.h \
    random.h \
    parameters.h \
    matingsystem.h \
    helperfunctions.h \
    femalesampling.h \
    densitydependentselection.h \
    enums.h
QT += core gui
TEMPLATE = app
RESOURCES +=
FORMS +=

 

 

 

 

 

ProjectChrisWiley.cpp

 

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
//---------------------------------------------------------------------------
USEFORM("UnitMain.cpp", FormMain);
USEFORM("UnitThreeDotsChasing.cpp", FormThreeDotsChasing);
USEFORM("UnitAboutBox2.cpp", FormAboutBox2);
//---------------------------------------------------------------------------
WINAPI WinMain(HINSTANCE, HINSTANCE, LPSTR, int)
{
  try
  {
     Application->Initialize();
     Application->Title = "The Chris Wiley Project";
                 Application->CreateForm(__classid(TFormMain), &FormMain);
     Application->Run();
  }
  catch (Exception &exception)
  {
     Application->ShowException(&exception);
  }
  catch (...)
  {
     try
     {
       throw Exception("");
     }
     catch (Exception &exception)
     {
       Application->ShowException(&exception);
     }
  }
  return 0;
}
//---------------------------------------------------------------------------

 

 

 

 

 

UnitAboutBox2.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef UnitAboutBox2H
#define UnitAboutBox2H
//---------------------------------------------------------------------------
#include <Classes.hpp>
#include <Controls.hpp>
#include <StdCtrls.hpp>
#include <Forms.hpp>
#include <ExtCtrls.hpp>
#include <Graphics.hpp>
//---------------------------------------------------------------------------
class TFormAboutBox2 : public TForm
{
__published: // IDE-managed Components
        TImage *Image1;
        TStaticText *StaticText1;
        TStaticText *StaticText2;
  TStaticText *StaticText3;
private: // User declarations
public: // User declarations
        __fastcall TFormAboutBox2(TComponent* Owner);
};
//---------------------------------------------------------------------------
extern PACKAGE TFormAboutBox2 *FormAboutBox2;
//---------------------------------------------------------------------------
#endif

 

 

 

 

 

UnitAboutBox2.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "UnitAboutBox2.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TFormAboutBox2 *FormAboutBox2;
//---------------------------------------------------------------------------
__fastcall TFormAboutBox2::TFormAboutBox2(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------

 

 

 

 

 

UnitMain.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef UnitMainH
#define UnitMainH
#include <Chart.hpp>
#include <Classes.hpp>
#include <ComCtrls.hpp>
#include <Controls.hpp>
#include <ExtCtrls.hpp>
#include <Grids.hpp>
#include <Series.hpp>
#include <StdCtrls.hpp>
#include <TeEngine.hpp>
#include <TeeProcs.hpp>
#include <AppEvnts.hpp>
//---------------------------------------------------------------------------
#include <mmsystem.h>
#include <Buttons.hpp>
#include <Dialogs.hpp>
//#include <Buttons.hpp>
//#include <utility>
//#include <functional>
//#include <iterator>
#include <string>
#include <vector>
#include <memory>
#include "UnitSimulation.h"
#include "UnitTimeSeries.h"
#include "UnitParameters.h"
//---------------------------------------------------------------------------
class TFormMain : public TForm
{
__published: // IDE-managed Components
  TStringGrid *StringGridParameters;
  TButton *ButtonRun;
  TStatusBar *StatusBarMain;
        TChart *ChartBias;
        TLineSeries *Series1;
        TLineSeries *Series2;
        TLineSeries *Series3;
        TPanel *PanelLeft;
        TRadioGroup *RadioGroupGamy;
        TRadioGroup *RadioGroupFemaleSampling;
        TProgressBar *ProgressBarSim;
  TPageControl *PageControlMain;
        TTabSheet *TabSheet1;
        TTabSheet *TabSheet3;
        TTabSheet *TabSheet4;
        TRichEdit *RichEdit2;
        TRadioGroup *RadioGroupDensityDependentSelection;
        TTabSheet *TabSheet7;
        TChart *ChartPopSize;
        TLineSeries *Series12;
        TLineSeries *Series13;
        TPanel *PanelLeftTop;
        TLineSeries *Series16;
        TLineSeries *Series19;
        TEdit *EditFractionFrom;
        TEdit *EditFractionTo;
        TEdit *EditFractionStep;
        TLabel *Label3;
        TLabel *Label4;
        TTabSheet *TabSheet6;
        TChart *ChartMateTime;
        TLineSeries *Series20;
        TLineSeries *Series21;
        TLineSeries *Series22;
        TLineSeries *Series23;
        TTabSheet *TabSheet2;
        TPanel *PanelTest;
        TChart *ChartTest;
        TStringGrid *StringGridTest;
        TButton *ButtonTest;
        TBarSeries *Series4;
        TTabSheet *TabSheet5;
        TChart *ChartTraitPreference;
        TPointSeries *Series5;
        TPointSeries *Series6;
        TLineSeries *Series7;
        TTabSheet *TabSheet12;
        TChart *ChartMateFraction;
        TLineSeries *Series9;
        TLineSeries *Series10;
        TLineSeries *Series11;
        TLineSeries *Series24;
        TPointSeries *Series25;
        TPanel *Panel1;
        TPageControl *PageControlSim;
        TTabSheet *TabSheet13;
        TTabSheet *TabSheet14;
        TLabel *Label1;
        TEdit *EditFraction;
        TLabel *Label2;
        TPageControl *PageControl2;
        TTabSheet *TabSheet8;
        TTabSheet *TabSheet9;
        TPanel *PanelTestProbTop;
        TButton *ButtonTestSurvivalSpecies;
        TChart *ChartTestProbabilities;
        TLineSeries *Series26;
        TBitBtn *BitBtn1;
        TTabSheet *TabSheet10;
        TChart *ChartBiasTime;
        TLineSeries *Series27;
        TLineSeries *Series28;
        TLineSeries *Series29;
        TButton *ButtonSave;
        TSaveDialog *SaveDialog1;
        TRichEdit *RichEditOutput;
        TButton *ButtonTestTrait;
        TButton *ButtonTestPreference;
        TButton *ButtonTestProbabilisticMating;
  TButton *ButtonTestHelp;
        TTabSheet *TabSheetAbout;
        TRichEdit *RichEditAbout;
  void __fastcall ButtonRunClick(TObject *Sender);
        void __fastcall ButtonTestClick(TObject *Sender);
        void __fastcall ButtonTestSurvivalSpeciesClick(TObject *Sender);
        void __fastcall FormMouseMove(TObject *Sender, TShiftState Shift,
          int X, int Y);
        void __fastcall BitBtn1Click(TObject *Sender);
        void __fastcall ButtonSaveClick(TObject *Sender);
        void __fastcall ButtonTestTraitClick(TObject *Sender);
        void __fastcall ButtonTestPreferenceClick(TObject *Sender);
        void __fastcall ButtonTestProbabilisticMatingClick(TObject *Sender);
  void __fastcall StringGridParametersSelectCell(TObject *Sender, int ACol,
          int ARow, bool &CanSelect);
  void __fastcall ButtonTestHelpClick(TObject *Sender);
  void __fastcall RadioGroupFemaleSamplingClick(TObject *Sender);
  void __fastcall RadioGroupDensityDependentSelectionClick(
          TObject *Sender);
  void __fastcall RadioGroupGamyClick(TObject *Sender);
  void __fastcall PageControlSimChange(TObject *Sender);
  void __fastcall BitBtn1MouseDown(TObject *Sender, TMouseButton Button,
          TShiftState Shift, int X, int Y);
private: // User declarations
  void setCursor(const String& cursorName);
  bool checkInput();
  void TFormMain::plotTimePoint(
    const TimePoint& mean,
    const TimePoint& stdError,
    const double& chartX);
  void TFormMain::plotTimeSeries(
    const TimeSeries& mean,
    const TimeSeries& stdError);
  Parameters readStringGrid() const;
public: // User declarations
  __fastcall TFormMain(TComponent* Owner);
};
//---------------------------------------------------------------------------
extern PACKAGE TFormMain *FormMain;
//---------------------------------------------------------------------------
std::string toString(const String& ansi);
String toAnsiString(const std::string& myString);
void emptyChart(TChart* chart);
std::auto_ptr<TStringList> getStringGrid(const TStringGrid * stringGrid);
std::auto_ptr<TStringList> getChart(const TChart * anyChart);
std::auto_ptr<TStringList> getChartSeries(const TChartSeries * series);


#endif

 

 

 

 

 

UnitMain.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "UnitMain.h"
#include "UnitAboutBox2.h"
#include "UnitThreeDotsChasing.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TFormMain *FormMain;
//---------------------------------------------------------------------------
__fastcall TFormMain::TFormMain(TComponent* Owner)
  : TForm(Owner)
{
  #ifdef NDEBUG
    OutputDebugString("NDEBUG defined. No output anymore");
    this->Caption = this->Caption + " NoDebug";
  #else
    OutputDebugString("NDEBUG not defined.");
    this->Caption = this->Caption + " Debug";
  #endif

  setCursor("CursorPiedFlycatcher.cur");
  StringGridTest->Cells[0][0] = "";
  StringGridTest->Cells[1][0] = "Female" ;
  StringGridTest->Cells[2][0] = "Male #1";
  StringGridTest->Cells[3][0] = "Male #2";
  StringGridTest->Cells[4][0] = "Male #3";
  StringGridTest->Cells[5][0] = "Male #4";
  StringGridTest->Cells[0][1] = "SpeciesValue";
  StringGridTest->Cells[1][1] = "-1.0" ;
  StringGridTest->Cells[2][1] = "-1.0";
  StringGridTest->Cells[3][1] = "-0.5";
  StringGridTest->Cells[4][1] = "0.5";
  StringGridTest->Cells[5][1] = "1.0";
  StringGridTest->Cells[0][2] = "Trait/Preference";
  StringGridTest->Cells[1][2] = "1.0" ;
  StringGridTest->Cells[2][2] = "-1.0";
  StringGridTest->Cells[3][2] = "-0.5";
  StringGridTest->Cells[4][2] = "0.5";
  StringGridTest->Cells[5][2] = "1.0";

  //Fill in StringGridParameters
  StringGridParameters->RowCount = 23;
  StringGridParameters->Cells[0][ 0] = "Parameter";             StringGridParameters->Cells[1][ 0] = "Value";
  StringGridParameters->Cells[0][ 1] = "Number of females";     StringGridParameters->Cells[1][ 1] = "1000";
  StringGridParameters->Cells[0][ 2] = "Number of males";       StringGridParameters->Cells[1][ 2] = "1000";
  StringGridParameters->Cells[0][ 3] = "Best-of-how-much";      StringGridParameters->Cells[1][ 3] = "5";
  StringGridParameters->Cells[0][ 4] = "Assessing error A";     StringGridParameters->Cells[1][ 4] = "0.001";
  StringGridParameters->Cells[0][ 5] = "Assessing error B";     StringGridParameters->Cells[1][ 5] = "0.001";
  StringGridParameters->Cells[0][ 6] = "Mean trait A";          StringGridParameters->Cells[1][ 6] = "-1.0";
  StringGridParameters->Cells[0][ 7] = "StdDev trait A";        StringGridParameters->Cells[1][ 7] = "0.0";
  StringGridParameters->Cells[0][ 8] = "Mean trait B";          StringGridParameters->Cells[1][ 8] = "1.0";
  StringGridParameters->Cells[0][ 9] = "StdDev trait B";        StringGridParameters->Cells[1][ 9] = "0.0";
  StringGridParameters->Cells[0][10] = "Mean preference A";     StringGridParameters->Cells[1][10] = "-1.0";
  StringGridParameters->Cells[0][11] = "StdDev preference A";   StringGridParameters->Cells[1][11] = "0.0";
  StringGridParameters->Cells[0][12] = "Mean preference B";     StringGridParameters->Cells[1][12] = "1.0";
  StringGridParameters->Cells[0][13] = "StdDev preference B";   StringGridParameters->Cells[1][13] = "0.0";
  StringGridParameters->Cells[0][14] = "Number of simulations"; StringGridParameters->Cells[1][14] = "1";
  StringGridParameters->Cells[0][15] = "Number of generations"; StringGridParameters->Cells[1][15] = "1";
  StringGridParameters->Cells[0][16] = "Number of offspring";   StringGridParameters->Cells[1][16] = "6";
  StringGridParameters->Cells[0][17] = "SurviveSpeciesAlpha";   StringGridParameters->Cells[1][17] = "1.0";
  StringGridParameters->Cells[0][18] = "SurviveSpeciesBeta";    StringGridParameters->Cells[1][18] = "10.0";
  StringGridParameters->Cells[0][19] = "SigmaSquared";          StringGridParameters->Cells[1][19] = "1.0";
  StringGridParameters->Cells[0][20] = "Cost trait";            StringGridParameters->Cells[1][20] = "0.1";
  StringGridParameters->Cells[0][21] = "Cost preference";       StringGridParameters->Cells[1][21] = "0.01";
  StringGridParameters->Cells[0][22] = "Mutation rate";         StringGridParameters->Cells[1][22] = "0.01";

  PageControlSimChange(0);
}
//---------------------------------------------------------------------------
void __fastcall TFormMain::ButtonRunClick(TObject *Sender)
{
  //Check the input in own scope
  if (checkInput()==false) return;

  StatusBarMain->Panels->Items[0]->Text = "Running simulation";
  setCursor("CursorPiedFlycatcherSleeping.cur");
  Refresh(); //To draw this text to the StatusBar

  emptyChart(ChartBias);
  emptyChart(ChartMateFraction);

  //Read parameters
  Parameters parameters = readStringGrid();
  parameters.matingSystem   = (RadioGroupGamy->ItemIndex  == 0 ? monogamy : polygyny);

  switch(RadioGroupFemaleSampling->ItemIndex)
  {
    case 0: parameters.femaleSampling = bestOfNconspicific;           break;
    case 1: parameters.femaleSampling = bestOfNextremeTrait;          break;
    case 2: parameters.femaleSampling = bestOfNclosestTrait;          break;
    case 3: parameters.femaleSampling = fixedThresholdConspicific;    break;
    case 4: parameters.femaleSampling = fixedThresholdTraitSign;      break;
    case 5: parameters.femaleSampling = fixedThresholdProbabilistic; break;
    default: assert(!"Unknown index of RadioGroupTestSampling"); std::exit(1);
  }

  parameters.densityDependentSelection = (RadioGroupDensityDependentSelection->ItemIndex == 0 ? afterMating : beforeMating );

  //Show parameters
  RichEditOutput->Lines->Clear();
  RichEditOutput->Lines->Add("Number of males: " + IntToStr(parameters.nMales));
  RichEditOutput->Lines->Add("Number of males: " + IntToStr(parameters.nMales));
  RichEditOutput->Lines->Add("Number of females: " + IntToStr(parameters.nFemales));
  RichEditOutput->Lines->Add("Best of how much: " + IntToStr(parameters.bestOfHowMuch));
  RichEditOutput->Lines->Add("Assessing error A: " + FloatToStr(parameters.assessingErrorA));
  RichEditOutput->Lines->Add("Assessing error B: " + FloatToStr(parameters.assessingErrorB));
  RichEditOutput->Lines->Add("Mean trait A: " + FloatToStr(parameters.meanTraitA));
  RichEditOutput->Lines->Add("StdDev trait A: " + FloatToStr(parameters.stdDevTraitA));
  RichEditOutput->Lines->Add("Mean trait males B: " + FloatToStr(parameters.meanTraitB));
  RichEditOutput->Lines->Add("StdDev trait males B: " + FloatToStr(parameters.stdDevTraitB));
  RichEditOutput->Lines->Add("Mean preference A: " + FloatToStr(parameters.meanPreferenceA));
  RichEditOutput->Lines->Add("StdDev preference A: " + FloatToStr(parameters.stdDevPreferenceA));
  RichEditOutput->Lines->Add("Mean preference males B: " + FloatToStr(parameters.meanPreferenceB));
  RichEditOutput->Lines->Add("StdDev preference males B: " + FloatToStr(parameters.stdDevPreferenceB));

  RichEditOutput->Lines->Add("Number of simulation: " + IntToStr(parameters.nSimulations));
  RichEditOutput->Lines->Add("Number of generations: " + IntToStr(parameters.nGenerations));
  RichEditOutput->Lines->Add("Number of offspring per couple: " + IntToStr(parameters.nOffspring));
  RichEditOutput->Lines->Add("SurviveSpeciesAlpha: " + FloatToStr(parameters.surviveSpeciesAlpha));
  RichEditOutput->Lines->Add("SurviveSpeciesBeta: " + FloatToStr(parameters.surviveSpeciesBeta));
  RichEditOutput->Lines->Add("SigmaSquared: " + FloatToStr(parameters.sigmaSquared));
  RichEditOutput->Lines->Add("Cost trait: " + FloatToStr(parameters.costTrait));
  RichEditOutput->Lines->Add("Cost preference: " + FloatToStr(parameters.costPreference));
  RichEditOutput->Lines->Add("Mutation rate: " + FloatToStr(parameters.mutationRate));
  RichEditOutput->Lines->Add(parameters.matingSystem  == monogamy ? "Mating system: monogamy" : "Mating system: polygyny" );
  switch(parameters.femaleSampling)
  {
    case bestOfNconspicific:
      RichEditOutput->Lines->Add("Female sampling: Best-Of-N conspicific");
      break;
    case bestOfNextremeTrait:
      RichEditOutput->Lines->Add("Female sampling: Best-Of-N extreme trait");
      break;
    case bestOfNclosestTrait:
      RichEditOutput->Lines->Add("Female sampling: Best-Of-N closest trait");
      break;
    case fixedThresholdConspicific:
      RichEditOutput->Lines->Add("Female sampling: Fixed threshold conspicific");
      break;
    case fixedThresholdTraitSign:
      RichEditOutput->Lines->Add("Female sampling: Fixed threshold trait sign");
      break;
    case fixedThresholdProbabilistic:
      RichEditOutput->Lines->Add("Female sampling: Fixed threshold probabilistic");
      break;
  }
  RichEditOutput->Lines->Add(parameters.densityDependentSelection==beforeMating ? "Density dependent selection: before mating" : "Density dependent selection: after mating");

  //Start simulation
  const double deltaFraction
    = (PageControlSim->ActivePageIndex == 1 ? EditFractionStep->Text.ToDouble() : 1.0);
  const double fractionFrom
    = (PageControlSim->ActivePageIndex == 1 ? EditFractionFrom->Text.ToDouble() : EditFraction->Text.ToDouble());
  const double fractionTo
    = (PageControlSim->ActivePageIndex == 1 ? EditFractionTo->Text.ToDouble() : 1.0);

  std::string lastErrorMessage;
  for (double fraction = fractionFrom; fraction < fractionTo; fraction+=deltaFraction)
  //double fraction = 0.9;
  {
    ProgressBarSim->Position = fraction * ProgressBarSim->Max;
    ProgressBarSim->Refresh();
    emptyChart(ChartBiasTime);
    emptyChart(ChartMateTime);
    emptyChart(ChartPopSize);
    emptyChart(ChartTraitPreference);

    parameters.fractionMaleA   = fraction;
    parameters.fractionFemaleA = fraction;
    //parameters.nGenerations = 1;
    const unsigned int nSimulations = parameters.nSimulations;
    std::vector<TimeSeries> allTimeSeries; //Size 0, use push_back
    for (unsigned int simulation = 0; simulation < nSimulations; ++simulation)
    {
      //Stick parameter in simulation and run it
      std::auto_ptr<Simulation> simulation(new Simulation(parameters));
      //simulation->showPopulation(StringGridDebug);
      simulation->execute();
      //Get the results from the sim
      //const TimeSeries thisSimTimeSeries = simulation->getTimeSeries();
      //timeSeries.push_back(thisSimTimeSeries);
      allTimeSeries.push_back(simulation->getTimeSeries());
      lastErrorMessage = simulation->mError;
    }
    //Plot the average of the results
    TimePoint timePointMean, timePointStdError;
    getMeanAndStdErrorEndPoint(allTimeSeries,timePointMean,timePointStdError);
    TimeSeries timeSeriesMean, timeSeriesStdError;
    getMeanAndStdErrorTimeSeries(allTimeSeries,timeSeriesMean, timeSeriesStdError);
    plotTimePoint(timePointMean, timePointStdError,fraction);
    plotTimeSeries(timeSeriesMean, timeSeriesStdError);
  }

  setCursor("CursorPiedFlycatcher.cur");
  StatusBarMain->Panels->Items[0]->Text = toAnsiString(lastErrorMessage);
  if (BitBtn1->Tag==0) PlaySound("PiedFlycatcher.wav",0,SND_FILENAME | SND_ASYNC);

}
//---------------------------------------------------------------------------

//Checks the input, returns true if all valid
bool TFormMain::checkInput()
{
  int checkInt; bool check; double checkDouble;
  //nFemales
  check = TryStrToInt(StringGridParameters->Cells[1][1],checkInt);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'nFemales' is not a valid input"; return false; }
  if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'nFemales' should be positive"; return false; }
  //nMales
  check = TryStrToInt(StringGridParameters->Cells[1][2],checkInt);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'nMales' is not a valid input"; return false; }
  if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'nMales' should be positive"; return false; }
  //bestOfHowMuch
  if (   RadioGroupFemaleSampling->ItemIndex == 0
      || RadioGroupFemaleSampling->ItemIndex == 1
      || RadioGroupFemaleSampling->ItemIndex == 2 ) //If female uses bestOfHowMuch
  {
    check = TryStrToInt(StringGridParameters->Cells[1][3],checkInt);
    if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'bestOfHowMuch' is not a valid input"; return false; }
    if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'bestOfHowMuch' should be bigger then 0"; return false; }
  }
  //assessingErrorA
  check = TryStrToFloat(StringGridParameters->Cells[1][4],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorA' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorA' should be positive"; return false; }
  if (checkDouble > 1.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorA' should be smaller or equal to 1.0"; return false; }
  //assessingErrorB
  check = TryStrToFloat(StringGridParameters->Cells[1][5],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorB' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorB' should be positive"; return false; }
  if (checkDouble > 1.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'assessingErrorB' should be smaller or equal to 1.0"; return false; }
  //Mean trait males A
  check = TryStrToFloat(StringGridParameters->Cells[1][6],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Mean trait males A' is not a valid input"; return false; }
  //StdDev trait males A
  check = TryStrToFloat(StringGridParameters->Cells[1][7],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'StdDev trait males A' is not a valid input"; return false; }
  //Mean trait males B
  check = TryStrToFloat(StringGridParameters->Cells[1][8],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Mean trait males B' is not a valid input"; return false; }
  //StdDev trait males B
  check = TryStrToFloat(StringGridParameters->Cells[1][9],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'StdDev trait males B' is not a valid input"; return false; }
  //Mean preference males A
  check = TryStrToFloat(StringGridParameters->Cells[1][10],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Mean preference males A' is not a valid input"; return false; }
  //StdDev preference males A
  check = TryStrToFloat(StringGridParameters->Cells[1][11],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'StdDev preference males A' is not a valid input"; return false; }
  //Mean preference males B
  check = TryStrToFloat(StringGridParameters->Cells[1][12],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Mean preference males B' is not a valid input"; return false; }
  //StdDev preference males B
  check = TryStrToFloat(StringGridParameters->Cells[1][13],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'StdDev preference males B' is not a valid input"; return false; }
  //Number of simulations
  check = TryStrToInt(StringGridParameters->Cells[1][14],checkInt);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of simulations' is not a valid input"; return false; }
  if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of simulations' should be positive"; return false; }
  //Number of generations
  check = TryStrToInt(StringGridParameters->Cells[1][15],checkInt);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of generations' is not a valid input"; return false; }
  if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of generations' should be positive"; return false; }
  //Number of offspring
  check = TryStrToInt(StringGridParameters->Cells[1][16],checkInt);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of offspring' is not a valid input"; return false; }
  if (checkInt <= 0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'Number of offspring' should be positive"; return false; }
  //SurviveSpeciesAlpha
  check = TryStrToFloat(StringGridParameters->Cells[1][17],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SurviveSpeciesAlpha' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SurviveSpeciesAlpha' should be positive"; return false; }
  if (checkDouble > 1.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SurviveSpeciesAlpha' should be smaller or equal to 1.0"; return false; }
  //SurviveSpeciesBeta
  check = TryStrToFloat(StringGridParameters->Cells[1][18],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SurviveSpeciesBeta' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SurviveSpeciesBeta' should be positive"; return false; }
  //SigmaSquared
  check = TryStrToFloat(StringGridParameters->Cells[1][19],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SigmaSquared' is not a valid input"; return false; }
  if (checkDouble <= 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'SigmaSquared' should be positive"; return false; }
  //CostTrait
  check = TryStrToFloat(StringGridParameters->Cells[1][20],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'CostTrait' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'CostTrait' should be positive"; return false; }
  //CostPreference
  check = TryStrToFloat(StringGridParameters->Cells[1][21],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'CostPreference' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'CostPreference' should be positive"; return false; }
  //MutationRate
  check = TryStrToFloat(StringGridParameters->Cells[1][22],checkDouble);
  if (check == false) { StatusBarMain->Panels->Items[0]->Text="Parameter 'MutationRate' is not a valid input"; return false; }
  if (checkDouble < 0.0) { StatusBarMain->Panels->Items[0]->Text="Parameter 'MutationRate' should be positive"; return false; }

  return true;
}
//---------------------------------------------------------------------------
void TFormMain::plotTimePoint(
  const TimePoint& mean,
  const TimePoint& stdError,
  const double& chartX)
{
  ChartBias->Series[0]->AddXY(chartX,mean.mateTally.calculateBiasA());
  ChartBias->Series[1]->AddXY(chartX,mean.mateTally.calculateBiasB());
  ChartBias->Series[2]->AddXY(chartX,mean.mateTally.calculateFractionMixedPairs());

  ChartMateFraction->Series[0]->AddXY(chartX,mean.mateTally.getNmateAA());
  ChartMateFraction->Series[1]->AddXY(chartX,mean.mateTally.getNmateAB());
  ChartMateFraction->Series[2]->AddXY(chartX,mean.mateTally.getNmateBA());
  ChartMateFraction->Series[3]->AddXY(chartX,mean.mateTally.getNmateBB());

}
//---------------------------------------------------------------------------
void TFormMain::plotTimeSeries(
  const TimeSeries& mean,
  const TimeSeries& stdError)
{
  emptyChart(ChartPopSize);
  emptyChart(ChartMateTime);
  const unsigned int nGenerations = mean.timePoints.size();
  for (unsigned i=0; i<nGenerations; ++i)
  {
    const double iD = static_cast<double>(i);
    /*
    ChartPopSize->Series[0]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNfemalesA());
    ChartPopSize->Series[1]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNfemaleHybridsA());
    ChartPopSize->Series[2]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNfemaleHybridsB());
    ChartPopSize->Series[3]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNfemalesB());
    ChartPopSize->Series[4]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNmalesA());
    ChartPopSize->Series[5]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNmaleHybridsA());
    ChartPopSize->Series[6]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNmaleHybridsB());
    ChartPopSize->Series[7]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNmalesB());
    ChartPopSize->Series[8]->AddXY(iD+0.00,mean.timePoints[i].speciesTallyOffspring.getNall());

    ChartPopSize->Series[0]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNfemalesA());
    ChartPopSize->Series[1]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNfemaleHybridsA());
    ChartPopSize->Series[2]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNfemaleHybridsB());
    ChartPopSize->Series[3]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNfemalesB());
    ChartPopSize->Series[4]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNmalesA());
    ChartPopSize->Series[5]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNmaleHybridsA());
    ChartPopSize->Series[6]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNmaleHybridsB());
    ChartPopSize->Series[7]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNmalesB());
    ChartPopSize->Series[8]->AddXY(iD+0.20,mean.timePoints[i].speciesTallyAfterSpeciesSelection.getNall());

    ChartPopSize->Series[0]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNfemalesA());
    ChartPopSize->Series[1]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNfemaleHybridsA());
    ChartPopSize->Series[2]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNfemaleHybridsB());
    ChartPopSize->Series[3]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNfemalesB());
    ChartPopSize->Series[4]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNmalesA());
    ChartPopSize->Series[5]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNmaleHybridsA());
    ChartPopSize->Series[6]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNmaleHybridsB());
    ChartPopSize->Series[7]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNmalesB());
    ChartPopSize->Series[8]->AddXY(iD+0.40,mean.timePoints[i].speciesTallyAfterTraitSelection.getNall());
    */
    ChartPopSize->Series[0]->AddXY(iD+0.60,mean.timePoints[i].speciesTallyAfterDensityDependentSelection.getNfemalesA());
    ChartPopSize->Series[1]->AddXY(iD+0.60,mean.timePoints[i].speciesTallyAfterDensityDependentSelection.getNfemalesB());
    ChartPopSize->Series[2]->AddXY(iD+0.60,mean.timePoints[i].speciesTallyAfterDensityDependentSelection.getNmalesA());
    ChartPopSize->Series[3]->AddXY(iD+0.60,mean.timePoints[i].speciesTallyAfterDensityDependentSelection.getNmalesB());
    //ChartPopSize->Series[4]->AddXY(iD+0.60,mean.timePoints[i].speciesTallyAfterDensityDependentSelection.getNall());

    ChartMateTime->Series[0]->AddXY(iD,mean.timePoints[i].mateTally.getNmateAA());
    ChartMateTime->Series[1]->AddXY(iD,mean.timePoints[i].mateTally.getNmateAB());
    ChartMateTime->Series[2]->AddXY(iD,mean.timePoints[i].mateTally.getNmateBA());
    ChartMateTime->Series[3]->AddXY(iD,mean.timePoints[i].mateTally.getNmateBB());
    ChartMateTime->Series[4]->AddXY(iD,mean.timePoints[i].mateTally.getNmateAll());

    ChartBiasTime->Series[0]->AddXY(iD,mean.timePoints[i].mateTally.calculateBiasA());
    ChartBiasTime->Series[1]->AddXY(iD,mean.timePoints[i].mateTally.calculateBiasB());
    ChartBiasTime->Series[2]->AddXY(iD,mean.timePoints[i].mateTally.calculateFractionMixedPairs());

    const unsigned int nTraits = mean.timePoints[i].traits.size();
    for (unsigned int j=0; j<nTraits; ++j) ChartTraitPreference->Series[0]->AddXY(iD,mean.timePoints[i].traits[j]);
    const unsigned int nPreferences = mean.timePoints[i].preferences.size();
    for (unsigned int j=0; j<nPreferences; ++j) ChartTraitPreference->Series[1]->AddXY(iD,mean.timePoints[i].preferences[j]);
    const unsigned int nSpeciesValues = mean.timePoints[i].descents.size();
    for (unsigned int j=0; j<nSpeciesValues; ++j) ChartTraitPreference->Series[2]->AddXY(iD,mean.timePoints[i].descents[j]);
  }
}
//---------------------------------------------------------------------------
void TFormMain::setCursor(const String& cursorName)
{
  TCursor myCursor = static_cast<TCursor>(22);
  Screen->Cursors[22] = LoadCursorFromFile(cursorName.c_str());
  Cursor = myCursor;
  FormMain->Cursor = myCursor;
  PageControlMain->Cursor = myCursor;
  StringGridParameters->Cursor = myCursor;
  PanelLeft->Cursor = myCursor;
  ChartBias->Cursor = myCursor;
  ChartPopSize->Cursor = myCursor;
  ChartMateTime->Cursor = myCursor;
  ButtonRun->Cursor = myCursor;
  ProgressBarSim->Cursor = myCursor;
  RadioGroupGamy->Cursor = myCursor;
  RadioGroupFemaleSampling->Cursor = myCursor;
  RichEditOutput->Cursor = myCursor;
  RichEdit2->Cursor = myCursor;
}
//---------------------------------------------------------------------------
Parameters TFormMain::readStringGrid() const
{
  Parameters parameters;
  parameters.nFemales            = StringGridParameters->Cells[1][ 1].ToInt();
  parameters.nMales              = StringGridParameters->Cells[1][ 2].ToInt();
  parameters.bestOfHowMuch       = StringGridParameters->Cells[1][ 3].ToInt();
  parameters.assessingErrorA     = StringGridParameters->Cells[1][ 4].ToDouble();
  parameters.assessingErrorB     = StringGridParameters->Cells[1][ 5].ToDouble();
  parameters.meanTraitA          = StringGridParameters->Cells[1][ 6].ToDouble();
  parameters.stdDevTraitA        = StringGridParameters->Cells[1][ 7].ToDouble();
  parameters.meanTraitB          = StringGridParameters->Cells[1][ 8].ToDouble();
  parameters.stdDevTraitB        = StringGridParameters->Cells[1][ 9].ToDouble();
  parameters.meanPreferenceA     = StringGridParameters->Cells[1][10].ToDouble();
  parameters.stdDevPreferenceA   = StringGridParameters->Cells[1][11].ToDouble();
  parameters.meanPreferenceB     = StringGridParameters->Cells[1][12].ToDouble();
  parameters.stdDevPreferenceB   = StringGridParameters->Cells[1][13].ToDouble();
  parameters.nSimulations        = StringGridParameters->Cells[1][14].ToInt();
  parameters.nGenerations        = StringGridParameters->Cells[1][15].ToInt();
  parameters.nOffspring          = StringGridParameters->Cells[1][16].ToInt();
  parameters.surviveSpeciesAlpha = StringGridParameters->Cells[1][17].ToDouble();
  parameters.surviveSpeciesBeta  = StringGridParameters->Cells[1][18].ToDouble();
  parameters.sigmaSquared        = StringGridParameters->Cells[1][19].ToDouble();
  parameters.costTrait           = StringGridParameters->Cells[1][20].ToDouble();
  parameters.costPreference      = StringGridParameters->Cells[1][21].ToDouble();
  parameters.mutationRate        = StringGridParameters->Cells[1][22].ToDouble();
  return parameters;
}
//---------------------------------------------------------------------------
void __fastcall TFormMain::ButtonTestClick(TObject *Sender)
{
  std::auto_ptr<StateFemaleSamplingBase> femaleSampling;
  switch(RadioGroupFemaleSampling->ItemIndex)
  {
    case 0: femaleSampling.reset(new StateFemaleSamplingBestOfNconspicific); break;
    case 1: femaleSampling.reset(new StateFemaleSamplingBestOfNextremeTrait); break;
    case 2: femaleSampling.reset(new StateFemaleSamplingBestOfNclosestTrait); break;
    case 3: femaleSampling.reset(new StateFemaleSamplingFixedThresholdConspicific); break;
    case 4: femaleSampling.reset(new StateFemaleSamplingFixedThresholdTraitSign); break;
    case 5: femaleSampling.reset(new StateFemaleSamplingFixedThresholdProbabilistic); break;
    default: assert(!"Unknown index of RadioGroupTestSampling"); std::exit(1);
  }

  Parameters parameters = readStringGrid();
  const double femaleSpeciesValue = StringGridTest->Cells[1][1].ToDouble();
  const double maleSpeciesValue1  = StringGridTest->Cells[2][1].ToDouble();
  const double maleSpeciesValue2  = StringGridTest->Cells[3][1].ToDouble();
  const double maleSpeciesValue3  = StringGridTest->Cells[4][1].ToDouble();
  const double maleSpeciesValue4  = StringGridTest->Cells[5][1].ToDouble();
  const double femalePreference   = StringGridTest->Cells[1][2].ToDouble();
  const double maleTrait1         = StringGridTest->Cells[2][2].ToDouble();
  const double maleTrait2         = StringGridTest->Cells[3][2].ToDouble();
  const double maleTrait3         = StringGridTest->Cells[4][2].ToDouble();
  const double maleTrait4         = StringGridTest->Cells[5][2].ToDouble();
  std::vector<Female> females = Bird::createTestFemales(parameters,femaleSpeciesValue, femalePreference);
  std::vector<Male> males = Bird::createTestMales(
    parameters,
    maleSpeciesValue1, maleTrait1,
    maleSpeciesValue2, maleTrait2,
    maleSpeciesValue3, maleTrait3,
    maleSpeciesValue4, maleTrait4
    );

  std::vector<unsigned int> histogramWinner(4);
  for (unsigned int i=0; i<4; ++i) histogramWinner[i]=0;

  for (unsigned int i=0; i<1000; ++i)
  {
    const unsigned int winnerIndex = femaleSampling->getWinnerIndex(males,females[0],parameters);
    ++(histogramWinner[winnerIndex]);
  }
  emptyChart(ChartTest);
  ChartTest->Series[0]->AddXY(0.0,histogramWinner[0]);
  ChartTest->Series[0]->AddXY(1.0,histogramWinner[1]);
  ChartTest->Series[0]->AddXY(2.0,histogramWinner[2]);
  ChartTest->Series[0]->AddXY(3.0,histogramWinner[3]);
}
//---------------------------------------------------------------------------
void __fastcall TFormMain::ButtonTestSurvivalSpeciesClick(TObject *Sender)
{
  ChartTestProbabilities->Series[0]->Clear();
  const double alpha = StringGridParameters->Cells[1][17].ToDouble();
  const double beta  = StringGridParameters->Cells[1][18].ToDouble();
  ChartTestProbabilities->Title->Text->Clear();
  ChartTestProbabilities->Title->Text->Add("const double y = 1.0 - (alpha * std::exp(-beta * descent * descent))");
  ChartTestProbabilities->Title->Text->Add("SurviveSpeciesAlpha: " + FloatToStr(alpha) + ", SurviveSpeciesBeta: " + FloatToStr(beta));
  ChartTestProbabilities->BottomAxis->Title->Caption = "SpeciesValue/descent";
  ChartTestProbabilities->LeftAxis->Title->Caption = "Chance of survival";
  for (double descent = -1.0; descent<1.0; descent+=0.01)
  {
    const double surviveSpecies = Simulation::chanceToSurviveSpecies(descent,alpha,beta);
    ChartTestProbabilities->Series[0]->AddXY(descent,surviveSpecies);
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::FormMouseMove(TObject *Sender,
      TShiftState Shift, int X, int Y)
{
  TCursor myCursor = static_cast<TCursor>(22);
  Screen->Cursors[22] = LoadCursorFromFile("CursorPiedFlycatcher.cur");
  Cursor = myCursor;
}
//---------------------------------------------------------------------------
void __fastcall TFormMain::ButtonSaveClick(TObject *Sender)
{
  std::auto_ptr<TStringList> stringList(new TStringList);

  stringList->Add(" ");
  stringList->Add("Parameters: ");

  if (PageControlSim->ActivePageIndex==0)
  {
    stringList->Add("FractionA : " + EditFraction->Text);
  }
  else
  {
    stringList->Add(
      "Fration A from: " + EditFractionFrom->Text
      + " to: " + EditFractionTo->Text
      + " in steps of: " + EditFractionStep->Text);
  }

  stringList->AddStrings(getStringGrid(StringGridParameters).release());

  switch(RadioGroupGamy->ItemIndex)
  {
    case 0: //monogamous
      stringList->Add("Mating system: monogamous");
      break;
    case 1: //polygynous
      stringList->Add("Mating system: polygynous");
      break;
  }

  switch(RadioGroupFemaleSampling->ItemIndex)
  {
    case 0: //bestOfNconspicific:
      stringList->Add("Female sampling: Best-Of-N conspicific");
      break;
    case 1: //bestOfNextremeTrait:
      stringList->Add("Female sampling: Best-Of-N extreme trait");
      break;
    case 2: //bestOfNclosestTrait:
      stringList->Add("Female sampling: Best-Of-N closest trait");
      break;
    case 3: //fixedThresholdConspicific:
      stringList->Add("Female sampling: Fixed threshold conspicific");
      break;
    case 4: //fixedThresholdTraitSign:
      stringList->Add("Female sampling: Fixed threshold trait sign");
      break;
    case 5: //fixedThresholdProbabilistic:
      stringList->Add("Female sampling: Fixed threshold probabilistic");
      break;
  }

  switch(RadioGroupDensityDependentSelection->ItemIndex)
  {
    case 0:
      stringList->Add("Density dependent selection: in reproduction, after mating");
      break;
    case 1:
      stringList->Add("Density dependent selection: after selection, before mating");
      break;
  }

  stringList->Add(" ");
  stringList->Add("Biases in time: ");
  stringList->Add("Time,BiasA,BiasB,Fraction of mixed pairs");
  stringList->AddStrings(getChart(ChartBiasTime).release());

  Dot("Trait and preferences in time");
  stringList->Add(" ");
  stringList->Add("Traits in time");
  stringList->Add("Time,Traits");
  stringList->AddStrings(getChartSeries(ChartTraitPreference->Series[0]).release());

  stringList->Add(" ");
  stringList->Add("Traits in time");
  stringList->Add("Time,Preference");
  stringList->AddStrings(getChartSeries(ChartTraitPreference->Series[1]).release());

  stringList->Add(" ");
  stringList->Add("Descent in time");
  stringList->Add("Time,Descent");
  stringList->AddStrings(getChartSeries(ChartTraitPreference->Series[2]).release());

  Dot("Number of matings in time");
  stringList->Add(" ");
  stringList->Add("Number of matings in time");
  stringList->Add("Time,AA,AB,BA,BB,Sum");
  stringList->AddStrings(getChart(ChartMateTime).release());

  Dot("number of individuals in time");
  stringList->Add(" ");
  stringList->Add("Number of individuals in time");
  stringList->Add("Time,FemalesA,FemalesB,MalesA,MalesB");
  stringList->AddStrings(getChart(ChartPopSize).release());

  Dot("Biases for different fractions of maleA: ");
  stringList->Add(" ");
  stringList->Add("Biases for different fractions of maleA: ");
  stringList->Add("Fraction,BiasA,BiasB,Fraction of mixed pairs");
  stringList->AddStrings(getChart(ChartBias).release());

  Dot("Number of matings for different proportions of maleA: ");
  stringList->Add(" ");
  stringList->Add("Number of matings for different proportions of maleA: ");
  stringList->Add("Fraction,AA,AB,BA,BB");
  stringList->AddStrings(getChart(ChartMateFraction).release());

  if (SaveDialog1->Execute()==true)
  {
    stringList->SaveToFile(SaveDialog1->FileName);
  }

  RichEditOutput->Lines = stringList.release();

}
//---------------------------------------------------------------------------
void __fastcall TFormMain::BitBtn1Click(TObject *Sender)
{
/*
  if (BitBtn1->Tag==1)
  {
    BitBtn1->Glyph->LoadFromFile("PiedFlycatcher.bmp");
    BitBtn1->Tag=0;
  }
  else
  {
    BitBtn1->Glyph->LoadFromFile("PiedFlycatcherSleeping.bmp");
    BitBtn1->Tag=1;
  }
  */
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
std::string toString(const String& ansi)
{
  const char * myChar = ansi.c_str();
  const std::string myString = myChar;
  return myString;
}
//---------------------------------------------------------------------------
String toAnsiString(const std::string& myString)
{
  const char * myChar = myString.c_str();
  const String myAnsi = myChar;
  return myAnsi;
}
//---------------------------------------------------------------------------
void emptyChart(TChart* chart)
{
  const int nSeries = chart->SeriesCount();
  for (int i=0; i<nSeries; ++i)
  {
    chart->Series[i]->Clear();
  }
}
//---------------------------------------------------------------------------
std::auto_ptr<TStringList> getStringGrid(const TStringGrid * stringGrid)
{
  //Those Borland people were not const correct, so I need a const cast here
  //Grumble, grumble...
  TStringGrid * grid = const_cast<TStringGrid*>(stringGrid);
  const String seperator = ",";

  std::auto_ptr<TStringList> stringList(new TStringList);
  const int maxy = grid->RowCount;
  const int maxx = grid->ColCount;
  for (int y=0; y<maxy; ++y)
  {
    String myString;
    for (int x=0; x<maxx-1; ++x)
    {
      myString+=grid->Cells[x][y] + seperator; //This function should have been const, as it is a read function
    }
    myString+=grid->Cells[maxx-1][y];
    stringList->Add(myString);
  }
  return stringList;
}
//---------------------------------------------------------------------------
std::auto_ptr<TStringList> getChart(const TChart * anyChart)
{
  //Those Borland people were not const correct, so I need a const cast here
  //Grumble, grumble...
  TChart * chart = const_cast<TChart*>(anyChart);
  const String seperator = ",";
  std::auto_ptr<TStringList> stringList(new TStringList);

  //Copy the pointers to the values
  std::vector< TChartValueList* > xValuesVector;
  std::vector< TChartValueList* > yValuesVector;
  const int nSeries = chart->SeriesCount();
  for (int i=0; i<nSeries; ++i)
  {
    xValuesVector.push_back(chart->Series[i]->XValues);
    yValuesVector.push_back(chart->Series[i]->YValues);
    //Assume there are as many X's as Y's
    assert(xValuesVector[i]->Count()==yValuesVector[i]->Count());
  }

  //Make a nice collumn of X values and then all Y values
  const int nRows = xValuesVector[0]->Count();
  #ifndef NDEBUG
    for (int i=0; i<nSeries; ++i)
    {
    assert(xValuesVector[i]->Count()==nRows);
    assert(yValuesVector[i]->Count()==nRows);
    }
  #endif

  const unsigned int nCols = xValuesVector.size();
  assert(nCols == yValuesVector.size());

  for (int y=0; y<nRows; ++y)
  {
    String myString = FloatToStr(xValuesVector[0]->operator [](y)) + seperator;
    for (unsigned int x=0; x<nCols-1; ++x)
    {
      myString+=FloatToStr(yValuesVector[x]->operator [](y)) + seperator;
    }
    myString+=FloatToStr(yValuesVector[nCols-1]->operator [](y));
    stringList->Add(myString);
  }

  return stringList;
}
//---------------------------------------------------------------------------
std::auto_ptr<TStringList> getChartSeries(const TChartSeries * series)
{
  std::auto_ptr<TStringList> stringList(new TStringList);
  const String seperator = ",";

  //Copy the pointers to the values
  TChartValueList * xValues = series->XValues;
  TChartValueList * yValues = series->YValues;
  //Assume there are as many X's as Y's
  assert(xValues->Count()==yValues->Count());

  const int nRows = xValues->Count();
  double x = -1.0;
  String myString;
  for (int y=0; y<nRows; ++y)
  {
    const double thisX = xValues->operator [](y);
    if (x!=thisX)
    {
      stringList->Add(myString);
      myString = FloatToStr(thisX) + seperator + FloatToStr(yValues->operator [](y)) + seperator;
      x = thisX;
    }
    else
    {
      myString += FloatToStr(yValues->operator [](y))+seperator;
    }
  }

  stringList->Add(myString);
  return stringList;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

void __fastcall TFormMain::ButtonTestTraitClick(TObject *Sender)
{
  ChartTestProbabilities->Series[0]->Clear();
  const double costTrait = StringGridParameters->Cells[1][20].ToDouble();
  ChartTestProbabilities->Title->Text->Clear();
  ChartTestProbabilities->Title->Text->Add("const double survival = std::exp(-costTrait * trait * trait);");
  ChartTestProbabilities->Title->Text->Add("costTrait: " + FloatToStr(costTrait));
  ChartTestProbabilities->BottomAxis->Title->Caption = "Species trait";
  ChartTestProbabilities->LeftAxis->Title->Caption = "Chance of survival";
  for (double maleTrait = -1.0; maleTrait<1.0; maleTrait+=0.01)
  {
    const double survivalTrait = Simulation::chanceToSurviveTrait(maleTrait,costTrait);
    ChartTestProbabilities->Series[0]->AddXY(maleTrait,survivalTrait);
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::ButtonTestPreferenceClick(TObject *Sender)
{
  ChartTestProbabilities->Series[0]->Clear();
  const double costPreference = StringGridParameters->Cells[1][21].ToDouble();

  ChartTestProbabilities->Title->Text->Clear();
  ChartTestProbabilities->Title->Text->Add("const double survival = std::exp(-costPreference * preference * preference");
  ChartTestProbabilities->Title->Text->Add("costPreference: " + FloatToStr(costPreference));
  ChartTestProbabilities->BottomAxis->Title->Caption = "Species preference";
  ChartTestProbabilities->LeftAxis->Title->Caption = "Chance of survival";
  for (double femalePreference = -1.0; femalePreference<1.0; femalePreference+=0.01)
  {
    const double survivalPreference = Simulation::chanceToSurvivePreference(femalePreference,costPreference);
    ChartTestProbabilities->Series[0]->AddXY(femalePreference,survivalPreference);
  }

}
//---------------------------------------------------------------------------

void __fastcall TFormMain::ButtonTestProbabilisticMatingClick(TObject *Sender)
{
  const double femalePreference = 0.0;
  const double sigmaSquared = StringGridParameters->Cells[1][19].ToDouble();
  ChartTestProbabilities->Series[0]->Clear();
  ChartTestProbabilities->Title->Text->Clear();
  ChartTestProbabilities->Title->Text->Add("const double chanceToMate =std::exp(-(maleTrait-femalePreference)*(maleTrait-femalePreference)/sigmaSquared);");
  ChartTestProbabilities->Title->Text->Add("FemalePreference: 0.0 (in this example)");
  ChartTestProbabilities->Title->Text->Add("sigmaSquared: " + FloatToStr(sigmaSquared));
  ChartTestProbabilities->BottomAxis->Title->Caption = "Male trait";
  ChartTestProbabilities->LeftAxis->Title->Caption = "Chance to mate";
  for (double maleTrait = -1.0; maleTrait<1.0; maleTrait+=0.01)
  {
    const double chanceToMate =
      StateFemaleSamplingFixedThresholdProbabilistic::getChanceToMate(femalePreference, maleTrait, sigmaSquared);
    ChartTestProbabilities->Series[0]->AddXY(maleTrait,chanceToMate);
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::StringGridParametersSelectCell(TObject *Sender,
      int ACol, int ARow, bool &CanSelect)
{
  //Show help
  switch (ARow)
  {
    case  1: StringGridParameters->Hint = "Number of females used in every generation"; break;
    case  2: StringGridParameters->Hint = "Number of males used in every generation"; break;
    case  3: StringGridParameters->Hint = "Number of males the female can sample from"; break;
    case  4: StringGridParameters->Hint = "Assessing error made by species A (0.0 = perfect sampling, 1.0=random mating)"; break;
    case  5: StringGridParameters->Hint = "Assessing error made by species B (0.0 = perfect sampling, 1.0=random mating)"; break;
    case  6: StringGridParameters->Hint = "Initial mean trait of species A"; break;
    case  7: StringGridParameters->Hint = "Initial StdDev of the trait of species A"; break;
    case  8: StringGridParameters->Hint = "Initial mean trait of species B"; break;
    case  9: StringGridParameters->Hint = "Initial StdDev of the trait of species B"; break;
    case 10: StringGridParameters->Hint = "Initial mean preference of species A"; break;
    case 11: StringGridParameters->Hint = "Initial StdDev of the preference of species A"; break;
    case 12: StringGridParameters->Hint = "Initial mean preference of species B"; break;
    case 13: StringGridParameters->Hint = "Initial StdDev of the preference of species B"; break;
    case 14: StringGridParameters->Hint = "Number of simulations of which the results are an average of"; break;
    case 15: StringGridParameters->Hint = "Number of generations"; break;
    case 16: StringGridParameters->Hint = "Number of offspring produced (used when density dependent selection = before reproduction)"; break;
    case 17: StringGridParameters->Hint = "Chance to die of pure hybrids (0.0 = Hybrids have equal survival, 1.0 = hybrids will die by chance 100%)"; break;
    case 18: StringGridParameters->Hint = "Increase in fitmess from hybrids to pure species (0.0 = hybrids and pure have equal survival, 100.0 = pure species have 100% survival)"; break;
    case 19: StringGridParameters->Hint = "Choosiness of female (only used in probabilistic mating, 1000.0 = random mating, 0.001 = very choosy"; break;
    case 20: StringGridParameters->Hint = "Cost trait in survival/trait (0.0 = trait is costless, 10.0 = trait is expensive)"; break;
    case 21: StringGridParameters->Hint = "Cost preference in survival/preference (0.0 = preference is costless, 10.0 = preference is expensive)"; break;
    case 22: StringGridParameters->Hint = "Standard deviation of mutation rate"; break;
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::ButtonTestHelpClick(TObject *Sender)
{
  ShowMessage(
    "These tests show the functions used."
    "\nThe results shown are derived from the functions used in the simulation"
    "\n(so no code duplication, except in the Chart titles)"
    "\nThe parameters used are from the Parameter StringGrid on the left."
    "\nNote that for these functions, the parameters are NOT checked."
    "\nTherefore, it is possible to give in a negative cost for trait,"
    "\nwhich will also yield a result, but when starting the simulation,"
    "\nthese values WILL be checked."
    );
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::RadioGroupFemaleSamplingClick(TObject *Sender)
{
  switch(RadioGroupFemaleSampling->ItemIndex)
  {
    case 0: //Best of N conspicific
      RadioGroupFemaleSampling->Hint = "The female chooses a conspicific from a lek of N males";
      break;
    case 1: //Best of N most extreme trait
      RadioGroupFemaleSampling->Hint = "The female chooses a male with the most extreme trait from a lek of N males. If her preference is smaller then 0.0, she prefers a male with trait smaller then 0.0";
      break;
    case 2: //Best of N closest trait-preference
      RadioGroupFemaleSampling->Hint = "The female chooses a male with trait closest near her preference from a lek of N males";
      break;
    case 3: //Fixed threshold conspicific
      RadioGroupFemaleSampling->Hint = "The female searches infinitely for a conspicific";
      break;
    case 4: //Fixed threshold same trait sign
      RadioGroupFemaleSampling->Hint = "The female searches infinitely for a male with the same trait sign as her preference (i.e. a negative-preference female prefers a negative-trait male)";
      break;
    case 5: //Fixed threshold probabilistic
      RadioGroupFemaleSampling->Hint = "The female searches infinitely for a male with trait close to her preference. The chance she will mate with him is defined by SigmaSquared";
      break;
  }
}
//---------------------------------------------------------------------------


void __fastcall TFormMain::RadioGroupDensityDependentSelectionClick(
      TObject *Sender)
{
  switch(RadioGroupDensityDependentSelection->ItemIndex)
  {
    case 0: //In reproduction, nOffspring = popSize parents
      RadioGroupDensityDependentSelection->Hint = "In reproduction, there are offspring produced until the original fractions of speciesA and speciesB are reached for a population size of nMales + nFemales";
    case 1: //After selection, nMatureOffspring = popSize parents
      RadioGroupDensityDependentSelection->Hint = "After selection, before mating, the population is cut down to the needed fractions of speciesA and speciesB to a population size of nMales + nFemales";
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::RadioGroupGamyClick(TObject *Sender)
{
  switch(RadioGroupGamy->ItemIndex)
  {
    case 0: //Monogamy
      RadioGroupGamy->Hint = "After a female has selected a male, they will be a couple for life (that is, one generation)";
      break;
    case 1: //Polygyny
      RadioGroupGamy->Hint = "A male can mate with multiple females";
      break;
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::PageControlSimChange(TObject *Sender)
{
  //ShowMessage("TEST" + IntToStr(PageControlSim->ActivePageIndex));
  if (PageControlSim->ActivePageIndex == 0)
  {
    //Single fraction sim (in time)
    PageControlMain->Pages[0]->TabVisible = true;
    PageControlMain->Pages[1]->TabVisible = true;
    PageControlMain->Pages[2]->TabVisible = true;
    PageControlMain->Pages[3]->TabVisible = true;
    PageControlMain->Pages[4]->TabVisible = false;
    PageControlMain->Pages[5]->TabVisible = false;
  }
  else
  {
    //Range of fractions sim
    PageControlMain->Pages[0]->TabVisible = false;
    PageControlMain->Pages[1]->TabVisible = false;
    PageControlMain->Pages[2]->TabVisible = false;
    PageControlMain->Pages[3]->TabVisible = false;
    PageControlMain->Pages[4]->TabVisible = true;
    PageControlMain->Pages[5]->TabVisible = true;
  }
}
//---------------------------------------------------------------------------

void __fastcall TFormMain::BitBtn1MouseDown(TObject *Sender,
      TMouseButton Button, TShiftState Shift, int X, int Y)
{
  if (BitBtn1->Tag==1)
  {
    BitBtn1->Glyph->LoadFromFile("PiedFlycatcher.bmp");
    BitBtn1->Tag=0;
  }
  else
  {
    BitBtn1->Glyph->LoadFromFile("PiedFlycatcherSleeping.bmp");
    BitBtn1->Tag=1;
  }
  if (Button == mbRight &&  X < 10 && Y < 10)
  {
    std::auto_ptr<TFormAboutBox2> aboutBox(new TFormAboutBox2(this));
    std::auto_ptr<TFormThreeDotsChasing> threeDotsChasing(new TFormThreeDotsChasing(this));
    //const int sumWidth = aboutBox->Width + threeDotsChasing->Width;
    aboutBox->Show();
    aboutBox->Left = 0;
    threeDotsChasing->ShowModal();
  }

}
//---------------------------------------------------------------------------

 

 

 

 

 

UnitThreeDotsChasing.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef UnitThreeDotsChasingH
#define UnitThreeDotsChasingH
//---------------------------------------------------------------------------
#include <Classes.hpp>
#include <Controls.hpp>
#include <StdCtrls.hpp>
#include <Forms.hpp>
#include <AppEvnts.hpp>
#include <ExtCtrls.hpp>
#include <Graphics.hpp>
//---------------------------------------------------------------------------
#define NDEBUG
#include <vector>
#include <assert>
#include <math>
#include <time>
//---------------------------------------------------------------------------
class TFormThreeDotsChasing : public TForm
{
__published: // IDE-managed Components
  TImage *Image1;
  TApplicationEvents *ApplicationEvents1;
  void __fastcall ApplicationEvents1Idle(TObject *Sender, bool &Done);
  void __fastcall FormClick(TObject *Sender);
  void __fastcall FormResize(TObject *Sender);
private: // User declarations
public: // User declarations
  __fastcall TFormThreeDotsChasing(TComponent* Owner);
  //int count;
  std::vector< std::vector<unsigned char> > distance;

  struct Point
  {
    short x;
    short y;
    char dx;
    char dy;
  };
  short abs(const short& i) { return (i<0 ? -i : i); }
  std::vector< Point > points;
  void resetPoints();
};
//---------------------------------------------------------------------------
extern PACKAGE TFormThreeDotsChasing *FormThreeDotsChasing;
//---------------------------------------------------------------------------
#endif

 

 

 

 

 

UnitThreeDotsChasing.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richèl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "UnitThreeDotsChasing.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TFormThreeDotsChasing *FormThreeDotsChasing;
//---------------------------------------------------------------------------
__fastcall TFormThreeDotsChasing::TFormThreeDotsChasing(TComponent* Owner)
  : TForm(Owner)
{
  points.resize(4);
  std::srand(clock());

  FormResize(0);
}
//---------------------------------------------------------------------------
void __fastcall TFormThreeDotsChasing::ApplicationEvents1Idle(TObject *Sender, bool &Done)
{
  //++count;
  //this->Caption = count;

  const char nPoints = points.size();
  const short maxy = Image1->Height;
  const short maxx = Image1->Width;
  for (char i=0; i<nPoints; ++i)
  {
    const char other = (i+1)%nPoints;
    assert(other>=0);
    assert(other<nPoints);
    if (points[i].x < points[other].x) ++(points[i].dx);
    else --(points[i].dx);
    if (points[i].y < points[other].y) ++(points[i].dy);
    else --(points[i].dy);
    points[i].x+=points[i].dx;
    if (points[i].x<0)
      { points[i].x=0; points[i].dx = -points[i].dx; }
    else if (points[i].x>=maxx)
      { points[i].x=maxx-1; points[i].dx = -points[i].dx; }
    points[i].y+=points[i].dy;
    if (points[i].y<0)
      { points[i].y=0; points[i].dy = -points[i].dy; }
    else if (points[i].y>=maxy)
      { points[i].y=maxy-1; points[i].dy = -points[i].dy; }
  }

  unsigned char * pLine;
  for (short y=0; y<maxy; ++y)
  {
    pLine=static_cast<unsigned char *>(Image1->Picture->Bitmap->ScanLine[y]);
    for (short x=0; x<maxx; ++x)
    {
      const short indexRedX   = abs(x-points[0].x)%maxx;
      const short indexRedY   = abs(y-points[0].y)%maxy;
      const short indexGreenX = abs(x-points[1].x)%maxx;
      const short indexGreenY = abs(y-points[1].y)%maxy;
      const short indexBlueX  = abs(x-points[2].x)%maxx;
      const short indexBlueY  = abs(y-points[2].y)%maxy;
      assert(indexRedX   >= 0 );
      assert(indexGreenX >= 0 );
      assert(indexBlueX  >= 0 );
      assert(indexRedY   >= 0 );
      assert(indexGreenY >= 0 );
      assert(indexBlueY  >= 0 );
      assert(indexRedX   < this->ClientWidth);
      assert(indexGreenX < this->ClientWidth);
      assert(indexBlueX  < this->ClientWidth);
      assert(indexRedY   < this->ClientHeight);
      assert(indexGreenY < this->ClientHeight);
      assert(indexBlueY  < this->ClientHeight);
      pLine[x*3+0]= distance[indexBlueX ][indexBlueY ]; //Blue
      pLine[x*3+1]= distance[indexGreenX][indexGreenY]; //Green
      pLine[x*3+2]= distance[indexRedX  ][indexRedY  ]; //Red
    }
  }
  this->Canvas->Draw(0,0,Image1->Picture->Graphic);
  Done = false;
}
//---------------------------------------------------------------------------
void __fastcall TFormThreeDotsChasing::FormClick(TObject *Sender)
{
  resetPoints();
}
//---------------------------------------------------------------------------
void TFormThreeDotsChasing::resetPoints()
{
  const short width  = this->ClientWidth;
  const short height = this->ClientHeight;
  const char nPoints = points.size();
  for (int i=0; i<nPoints; ++i)
  {
    points[i].x = std::rand()%width;
    points[i].y = std::rand()%height;
    points[i].dx = (std::rand()%5)-2;
    points[i].dy = (std::rand()%5)-2;
  }
}
//---------------------------------------------------------------------------
void __fastcall TFormThreeDotsChasing::FormResize(TObject *Sender)
{
  //Generate vectors for look-up table
  const short width  = this->ClientWidth;
  const short height = this->ClientHeight;

  distance.resize(width);
  for (int i=0; i<width; ++i) distance[i].resize(height);

  //Generate look-up table
  const double widthD  = static_cast<double>(width);
  const double heightD = static_cast<double>(height);
  const double maxDist = sqrt( (widthD*widthD) + (heightD*heightD) );
  for (int x=0; x<width; ++x)
  {
    for (int y=0; y<height; ++y)
    {
      const double xD = static_cast<double>(x);
      const double yD = static_cast<double>(y);
      const double dist = sqrt( (xD*xD) + (yD*yD) );
      //distance[x][y] = 256.0 * dist / maxDist;
      distance[x][y] = 255.0 * std::exp(-dist/150.0);
    }
  }
  //Points
  resetPoints();
  //Create image
  Image1->Picture->Bitmap->Width  = width;
  Image1->Picture->Bitmap->Height = height;

}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

bird.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef BIRD_H
#define BIRD_H
//---------------------------------------------------------------------------
#include <iostream>
#include <vector>
#include <cassert>
#include "helperfunctions.h"
#include "parameters.h"
#include "enums.h"
#include "random.h"

class Couple; //Forward declaration

class Bird
{
  public:

  Bird(const Bird& mother, const Bird& father, const double& mutation);

  static std::vector<Bird> createMales(const Parameters& parameters);
  static std::vector<Bird> createFemales(const Parameters& parameters);

  static std::vector<Bird> createTestMales(
    const Parameters& parameters,
    const double& species1, const double& trait1,
    const double& species2, const double& trait2,
    const double& species3, const double& trait3,
    const double& species4, const double& trait4);

  static std::vector<Bird> createTestFemales(
    const Parameters& parameters,
    const double& species, const double& preference);

  unsigned int index;
  double descent;
  enumSpecies species;
  double trait;
  double preference;
  double assessingError;


  private:
  //No empty constructor, Birds are either created from parents or
  //from the static createMales/createFemales functions
  Bird() {}
};

typedef Bird Male;
typedef Bird Female;

struct Offspring
{
  std::vector<Female> females;
  std::vector<Male>   males;
};

typedef Offspring Population;

class Couple
{
  public:
  Couple(const Female& anyFemale, const Male& anyMale)
    : female(anyFemale),
    male(anyMale)
  {
    //Nothing
  }
  Female female;
  Male   male;
};
//---------------------------------------------------------------------------
#endif // BIRD_H

 

 

 

 

 

bird.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "bird.h"
#include <cstdlib>
//---------------------------------------------------------------------------
Bird::Bird(const Bird& mother, const Bird& father, const double& mutation)
{
  descent = (mother.descent + father.descent) / 2.0;
  assert(descent>=-1.0 && descent<=1.0);
  const double chanceBeingPied = 1.0-((descent+1.0)/2.0);
  assert(chanceBeingPied>=0.0 && chanceBeingPied<=1.0);
  species = (rnd::uniform() > chanceBeingPied ? collaredFlycatcher : piedFlycatcher);
  //Make a bit
  unsigned int number = std::rand()%64;
  index = (number%2==0 ? mother.index : father.index);
  number>>=1;
  trait = (number%2==0 ? mother.trait : father.trait);
  trait+=rnd::normal(0.0,mutation);
  number>>=1;
  preference = (number%2==0 ? mother.preference : father.preference);
  preference+=rnd::normal(0.0,mutation);
  number>>=1;
  assessingError = (number%2==0 ? mother.assessingError : father.assessingError);

}
//---------------------------------------------------------------------------
std::vector<Male> Bird::createMales(const Parameters& parameters)
{
  //Create males
  assert(parameters.fractionMaleA>=0.0 && parameters.fractionMaleA<=1.0);
  //Create an empty vector of males
  std::vector<Male> males; //(parameters.nMales);
  //Calculate the number of males of speciesA
  const unsigned int nSpeciesA =
    0.5 + (static_cast<double>(parameters.nMales) * parameters.fractionMaleA)
    * (parameters.densityDependentSelection == beforeMating ? 2.0 : 1.0);
  const unsigned int nSpeciesB =
    0.5 + (static_cast<double>(parameters.nMales) * (1.0-parameters.fractionMaleA))
    * (parameters.densityDependentSelection == beforeMating ? 2.0 : 1.0);
  //Dot("Of " + IntToStr(nMales) + " males, there are " + IntToStr(nSpeciesA) + " of speciesA");
  //Fill the vector
  for (unsigned int i=0; i<nSpeciesA; ++i)
  {
    Male maleA;
    maleA.index = i;
    maleA.species = piedFlycatcher;
    maleA.descent = -1.0;
    maleA.trait = rnd::normal(parameters.meanTraitA, parameters.stdDevTraitA);
    maleA.preference = rnd::normal(parameters.meanPreferenceA, parameters.stdDevPreferenceA);
    maleA.assessingError = parameters.assessingErrorA;
    males.push_back(maleA);
  }
  for (unsigned int i=0 ; i<nSpeciesB; ++i)
  {
    Male maleB;
    maleB.index = i;
    maleB.species = collaredFlycatcher;
    maleB.descent = 1.0;
    maleB.trait = rnd::normal(parameters.meanTraitB, parameters.stdDevTraitB);
    maleB.preference = rnd::normal(parameters.meanPreferenceB, parameters.stdDevPreferenceB);
    maleB.assessingError = parameters.assessingErrorB;
    males.push_back(maleB);
  }
  return males;
}
//---------------------------------------------------------------------------
std::vector<Female> Bird::createFemales(const Parameters& parameters)
//Create females
{
  assert(parameters.fractionFemaleA>=0.0 && parameters.fractionFemaleA<=1.0);
  //Create an empty vector of females
  std::vector<Female> females;
  //Calculate the number of females of speciesA
  const unsigned int nSpeciesA
    = 0.5 + (static_cast<double>(parameters.nFemales) * parameters.fractionFemaleA)
    * (parameters.densityDependentSelection == beforeMating ? 2.0 : 1.0);

  const unsigned int nSpeciesB
    = 0.5 + (static_cast<double>(parameters.nFemales) * (1.0-parameters.fractionFemaleA))
    * (parameters.densityDependentSelection == beforeMating ? 2.0 : 1.0);

  //Fill the vector
  for (unsigned int i=0; i<nSpeciesA; ++i)
  {
    Female femaleA;
    femaleA.index = i;
    femaleA.species = piedFlycatcher;
    femaleA.descent = -1.0;
    femaleA.preference = rnd::normal(parameters.meanPreferenceA, parameters.stdDevPreferenceA);
    femaleA.assessingError = parameters.assessingErrorA;
    femaleA.trait = rnd::normal(parameters.meanTraitA, parameters.stdDevTraitA);
    females.push_back(femaleA);
  }
  for (unsigned int i=0 ; i<nSpeciesB; ++i)
  {
    Female femaleB;
    femaleB.index = i;
    femaleB.species = collaredFlycatcher;
    femaleB.descent = 1.0;
    femaleB.preference = rnd::normal(parameters.meanPreferenceB, parameters.stdDevPreferenceB);
    femaleB.assessingError = parameters.assessingErrorB;
    femaleB.trait = rnd::normal(parameters.meanTraitB, parameters.stdDevTraitB);
    females.push_back(femaleB);
  }
  return females;
}
//---------------------------------------------------------------------------
std::vector<Male> Bird::createTestMales(
  const Parameters& /*parameters*/,
  const double& species1, const double& trait1,
  const double& species2, const double& trait2,
  const double& species3, const double& trait3,
  const double& species4, const double& trait4
  )
{
  //Create an empty vector of males
  std::vector<Male> males;
  //Calculate the number of males of speciesA
  //Dot("Of " + IntToStr(nMales) + " males, there are " + IntToStr(nSpeciesA) + " of speciesA");
  //Fill the vector
  Male male1, male2, male3, male4;
  male1.descent = species1; male1.trait = trait1; male1.species = (rnd::uniform() > 1.0-((male1.descent+1.0)/2.0) ? collaredFlycatcher : piedFlycatcher);
  male2.descent = species2; male2.trait = trait2; male2.species = (rnd::uniform() > 1.0-((male2.descent+1.0)/2.0) ? collaredFlycatcher : piedFlycatcher);
  male3.descent = species3; male3.trait = trait3; male3.species = (rnd::uniform() > 1.0-((male3.descent+1.0)/2.0) ? collaredFlycatcher : piedFlycatcher);
  male4.descent = species4; male4.trait = trait4; male4.species = (rnd::uniform() > 1.0-((male4.descent+1.0)/2.0) ? collaredFlycatcher : piedFlycatcher);
  males.push_back(male1);
  males.push_back(male2);
  males.push_back(male3);
  males.push_back(male4);
  return males;
}
//---------------------------------------------------------------------------
std::vector<Female> Bird::createTestFemales(
    const Parameters& parameters,
    const double& species, const double& preference)
//Create females
{
  //Create an empty vector of females
  std::vector<Female> females;
  //Fill the vector with one female
  Female female;
  female.descent = species;
  female.species = (rnd::uniform() > 1.0-((female.descent+1.0)/2.0) ? collaredFlycatcher : piedFlycatcher);
  female.preference = preference;
  female.assessingError = parameters.assessingErrorA;
  females.push_back(female);
  return females;
}
//---------------------------------------------------------------------------

 

 

 

 

 

densitydependentselection.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef DENSITYDEPENDENTSELECTION_H
#define DENSITYDEPENDENTSELECTION_H
//---------------------------------------------------------------------------
#include <vector>
#include <string>
#include <algorithm>
//#include "UnitTallies.h"
#include "parameters.h"
#include "bird.h"

class StateDensityDependentSelectionBase
{
  public:
  StateDensityDependentSelectionBase() {}
  virtual std::string getString() const = 0;
  virtual bool canDoSelection(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    std::string& errorMessage
    //const SpeciesTally& speciesTally
    ) const = 0;
  virtual void doSelection(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters
    //const SpeciesTally& speciesTally
    ) const = 0;
};
//---------------------------------------------------------------------------
//class StateDensityDependentSelectionAfterMating : public StateDensityDependentSelectionBase
//{
//  public:
//  StateDensityDependentSelectionAfterMating() {}
//  std::string getString() const { return "DensityDependentSelectionAfterMating"; }
//  bool canDoSelection(
//    std::vector<Female>& females,
//    std::vector<Male>& males,
//    const Parameters& parameters,
//    std::string& errorMessage
//    //const SpeciesTally& speciesTally
//    ) const;
//  void doSelection(
//    std::vector<Female>& females,
//    std::vector<Male>& males,
//    const Parameters& parameters
//    //const SpeciesTally& speciesTally
//    ) const;
//};
//---------------------------------------------------------------------------
class StateDensityDependentSelectionBeforeMating : public StateDensityDependentSelectionBase
{
  public:
  StateDensityDependentSelectionBeforeMating() {}
  std::string getString() const { return "DensityDependentSelectionBeforeMating"; }
  bool canDoSelection(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    std::string& errorMessage
    //const SpeciesTally& speciesTally
    ) const;
  void doSelection(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters
    //const SpeciesTally& speciesTally
    ) const;
};
//---------------------------------------------------------------------------
#endif // DENSITYDEPENDENTSELECTION_H

 

 

 

 

 

densitydependentselection.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "densitydependentselection.h"
#include "speciestally.h"
//---------------------------------------------------------------------------
//bool StateDensityDependentSelectionAfterMating::canDoSelection(
//  std::vector<Female>& /*females*/,
//  std::vector<Male>& /*males*/,
//  const Parameters& /*parameters*/,
//  std::string& /*errorMessage*/
//  //const SpeciesTally& speciesTally
//  ) const
//{
//  return true;
//}
//---------------------------------------------------------------------------
//void StateDensityDependentSelectionAfterMating::doSelection(
//  std::vector<Female>& /*females*/,
//  std::vector<Male>& /*males*/,
//  const Parameters& /*parameters*/
//  //const SpeciesTally& speciesTally
//  ) const
//{
//  //Nothing, density dependent selection is found in MatingSystemFixedNumberOffspring
//}
//---------------------------------------------------------------------------
//This type of selection occurs before mating
//Idea is, that the population structure is constant before mating.
//This means that the number of malesA, malesB, femalesA and femalesB
//is constant. This is reached by removing superfluous individuals.
//If however, too few individuals are present, selection cannot take
//place and the population is called extinct.
bool StateDensityDependentSelectionBeforeMating::canDoSelection(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters,
  std::string& errorMessage
  //const SpeciesTally& speciesTally
  ) const
{
  //assert(speciesTally.isNull()==false);
  SpeciesTally speciesTally;
  speciesTally.tallySpecies(females,males);

  if ( speciesTally.getNfemalesA() < parameters.getNfemalesAwanted())
  {
    errorMessage = "Too few femalesA for selection."
      " Needed: " + itoa(parameters.getNfemalesAwanted())
      + " ,got: " + itoa(speciesTally.getNfemalesA());
    Dot(errorMessage);
    return false;
  }
  if ( speciesTally.getNfemalesB() < parameters.getNfemalesBwanted())
  {
    errorMessage = "Too few femalesB for selection."
      " Needed: " + itoa(parameters.getNfemalesBwanted())
      + " ,got: " + itoa(speciesTally.getNfemalesB());
    Dot(errorMessage);
    return false;
  }
  if ( speciesTally.getNmalesA() < parameters.getNmalesAwanted())
  {
    errorMessage = "Too few malesA for selection."
      " Needed: " + itoa(parameters.getNmalesAwanted())
      + " ,got: " + itoa(speciesTally.getNmalesA());
    Dot(errorMessage);
    return false;
  }
  if ( speciesTally.getNmalesB() < parameters.getNmalesBwanted())
  {
    errorMessage =  "Too few malesB for selection."
      " Needed: " + itoa(parameters.getNmalesBwanted())
      + " ,got: " + itoa(speciesTally.getNmalesB());
    Dot(errorMessage);
    return false;
  }

  return true;
}
//---------------------------------------------------------------------------
void StateDensityDependentSelectionBeforeMating::doSelection(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters
  //const SpeciesTally& speciesTally
  ) const
{
  unsigned int nMalesAwanted   = parameters.getNmalesAwanted();
  unsigned int nMalesBwanted   = parameters.getNmalesBwanted();
  assert(nMalesAwanted + nMalesBwanted == parameters.nMales);
  unsigned int nFemalesAwanted = parameters.getNfemalesAwanted();
  unsigned int nFemalesBwanted = parameters.getNfemalesBwanted();
  assert(nFemalesAwanted + nFemalesBwanted == parameters.nFemales);

  for (unsigned int i=0; i<males.size(); ) //No increment
  {
    if (males[i].species == piedFlycatcher) //Male is of type A
    {
      if (nMalesAwanted==0)
      {
        //Erase the male
         males.erase(males.begin() + i);
      }
      else
      {
        //Keep the male
        --nMalesAwanted;
        ++i;
      }
    }
    else //Male is of type B
    {
      if (nMalesBwanted==0)
      {
        //Erase the male
         males.erase(males.begin() + i);
      }
      else
      {
        //Keep the male
        --nMalesBwanted;
        ++i;
      }
    }
  }

  for (unsigned int i=0; i<females.size(); ) //No increment
  {
    if (females[i].species == piedFlycatcher) //Female is of type A
    {
      if (nFemalesAwanted==0)
      {
        //Erase the female
         females.erase(females.begin() + i);
      }
      else
      {
        //Keep the female
        --nFemalesAwanted;
        ++i;
      }
    }
    else //Female is of type B
    {
      if (nFemalesBwanted==0)
      {
        //Erase the female
         females.erase(females.begin() + i);
      }
      else
      {
        //Keep the female
        --nFemalesBwanted;
        ++i;
      }
    }
  }
  #ifdef __HACK_089236408766433745
    const unsigned int fs = females.size();
    const unsigned int ms = males.size();
    const unsigned int fw = parameters.nFemales;
    const unsigned int mw = parameters.nMales;
  #endif
  std::random_shuffle(females.begin(), females.end());
  std::random_shuffle(males.begin(), males.end());
  assert(females.size()==parameters.nFemales);
  assert(males.size()==parameters.nMales);
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

enums.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef ENUMS_H
#define ENUMS_H
//---------------------------------------------------------------------------
enum enumSpecies
{
  piedFlycatcher, //descent/speciesValue = -1
  collaredFlycatcher  //descent/speciesValue = 1.0
};

enum enumDensityDependentSelection
{
  beforeMating,
  afterMating
};

enum enumMatingSystem
{
  monogamy,
  polygyny //Male can mate multiple times
};

enum enumFemaleSampling
{
  bestOfNconspicific,
  bestOfNextremeTrait,
  bestOfNclosestTrait,
  fixedThresholdConspicific,
  fixedThresholdTraitSign,
  fixedThresholdProbabilistic
};


#endif // ENUMS_H

 

 

 

 

 

femalesampling.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef FEMALESAMPLING_H
#define FEMALESAMPLING_H
//---------------------------------------------------------------------------
#include <vector>
#include <string>
#include "bird.h"
#include "parameters.h"
#include "random.h"

//---------------------------------------------------------------------------
class StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingBase() {}
  virtual unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const = 0;
  virtual std::string getFemaleSampling() const = 0;
};
//---------------------------------------------------------------------------
class StateFemaleSamplingBestOfNconspicific : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingBestOfNconspicific() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Best of N conspicific"; }
};
//---------------------------------------------------------------------------
class StateFemaleSamplingBestOfNextremeTrait : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingBestOfNextremeTrait() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Best of N extreme trait"; }
};
//---------------------------------------------------------------------------
class StateFemaleSamplingBestOfNclosestTrait : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingBestOfNclosestTrait() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Best of N closest trait"; }
};
//---------------------------------------------------------------------------
class StateFemaleSamplingFixedThresholdConspicific : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingFixedThresholdConspicific() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Fixed threshold conspicific"; }
};
//---------------------------------------------------------------------------
class StateFemaleSamplingFixedThresholdTraitSign : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingFixedThresholdTraitSign() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Fixed threshold trait sign"; }
};
//---------------------------------------------------------------------------
class StateFemaleSamplingFixedThresholdProbabilistic : public StateFemaleSamplingBase
{
  public:
  StateFemaleSamplingFixedThresholdProbabilistic() {}
  unsigned int getWinnerIndex(const std::vector<Male>& males, const Female& female, const Parameters& parameters) const;
  std::string getFemaleSampling() const { return "Fixed threshold probabilistic"; }
  static inline double getChanceToMate(const double& femalePreference, const double& maleTrait, const double& sigmaSquared);


};
//---------------------------------------------------------------------------
std::vector<unsigned int> createIndices(const unsigned int& bestOfHowMuch, const unsigned int& nMales);
std::vector<unsigned int> createIndicesUnique(const unsigned int& bestOfHowMuch, const unsigned int& nMales);
//std::vector<double> createSpeciesVector(const std::vector<Male>& males, const std::vector<unsigned int>& maleIndex);
//std::vector<double> createTraitsVector(const std::vector<Male>& males, const std::vector<unsigned int>& maleIndex);
//unsigned int getBestMaleIndex(const std::vector<double>& maleTraits, const double& femalePreference, const double& assessingError);
//---------------------------------------------------------------------------
#endif // FEMALESAMPLING_H

 

 

 

 

 

femalesampling.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <algorithm>
#include "femalesampling.h"
//---------------------------------------------------------------------------
#define EXTREME_TEST
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingBestOfNconspicific::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& parameters) const
{
  const unsigned int bestOfHowMuch = parameters.bestOfHowMuch;
  //Create the male indices
  std::vector<unsigned int> maleIndices = createIndicesUnique(bestOfHowMuch, males.size());
  const std::vector<unsigned int> maleIndicesCopy = maleIndices;

  //Remove the males that are either heterospecific or assessed wrong
  for (unsigned int i=0; i<maleIndices.size(); ) //No ++i indeed
  {
    if (rnd::uniform() > female.assessingError) //The female does NOT make an error
    {
      if (female.species == piedFlycatcher) //Female is of species A (piedFlyCatcher)
      {
        if (males[ maleIndices[i] ].species == collaredFlycatcher) //Male is of species B
        {
          maleIndices.erase(maleIndices.begin() + i);
        }
        else //Male is of species A (piedFlycatcher)
          { ++i; } //The male is a conspicific and remains in the vecor
      }
      else //Female is of species B (collaredFlycatcher)
      {
        if (males[ maleIndices[i] ].species == collaredFlycatcher) //Male is of species B
          { ++i; }//The male is a conspicific and remains in the vecor
        else //Male is of species A
        {
          maleIndices.erase(maleIndices.begin() + i);
          //maleIndices.erase(&maleIndices.at(i));
        }
      }
    }
    else
    {
      ++i; //The wrongly assessed male remains
    }
  }

  //Dot("Remaining male Indices"); Dot(maleIndices);

  //With the remaining candidats, she mates at random
  if (maleIndices.size()>0) //She found a suitable candidat
  {
    const unsigned int candidatIndex = std::rand()%maleIndices.size();
    const unsigned int winnerIndex = maleIndices[candidatIndex];
    return winnerIndex;
  }
  else
  {
    return maleIndicesCopy[ std::rand()%maleIndicesCopy.size() ];
  }
}
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingBestOfNextremeTrait::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& parameters) const
{
  const unsigned int bestOfHowMuch = ( parameters.bestOfHowMuch < males.size() ? parameters.bestOfHowMuch : males.size());
  //Create the male indices
  std::vector<unsigned int> maleIndices = createIndicesUnique(bestOfHowMuch, males.size());

  //Dot("Male Indices"); Dot(maleIndices);

  if (rnd::uniform() > female.assessingError) //The female does NOT make an error
  { //Get the best male
    if (female.preference < 0.0) //Find male with lowest trait
    {
      unsigned int candidatIndex = 0;
      double bestTrait = males[ maleIndices[0] ].trait;
      for (unsigned int i=1; i<bestOfHowMuch; ++i)
      {
        if (males[ maleIndices[i] ].trait < bestTrait)
        {
          bestTrait = males[ maleIndices[i] ].trait;
          candidatIndex = i;
        }
      }
      const unsigned int winnerIndex = maleIndices[candidatIndex];
      return winnerIndex;
    }
    else //Find male with highest trait
    {
      unsigned int candidatIndex = 0;
      double bestTrait = males[ maleIndices[0] ].trait;
      for (unsigned int i=1; i<bestOfHowMuch; ++i)
      {
        if (males[ maleIndices[i] ].trait > bestTrait)
        {
          bestTrait = males[ maleIndices[i] ].trait;
          candidatIndex = i;
        }
      }
      const unsigned int winnerIndex = maleIndices[candidatIndex];
      return winnerIndex;
    }
  }
  else
  {
    //Female mates at random
    const unsigned int candidatIndex = std::rand()%bestOfHowMuch;
    unsigned int winnerIndex = maleIndices[candidatIndex];
    return winnerIndex;
  }
}
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingBestOfNclosestTrait::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& parameters) const
{
  const unsigned int bestOfHowMuch = ( parameters.bestOfHowMuch < males.size() ? parameters.bestOfHowMuch : males.size());
  //Create the male indices
  std::vector<unsigned int> maleIndices = createIndicesUnique(bestOfHowMuch, males.size());

  if (rnd::uniform() > female.assessingError) //The female does NOT make an error
  { //Get the best male
    unsigned int candidatIndex = 0;
    double closestTrait = std::fabs(males[ maleIndices[0] ].trait - female.preference);
    for (unsigned int i=1; i<bestOfHowMuch; ++i)
    {
      if (std::fabs(males[ maleIndices[i] ].trait - female.preference) < closestTrait)
      {
        closestTrait = std::fabs(males[ maleIndices[0] ].trait - female.preference);
        candidatIndex = i;
      }
    }
    const unsigned int winnerIndex = maleIndices[candidatIndex];
    return winnerIndex;
  }
  else
  {
    //Female mates at random
    const unsigned int candidatIndex = std::rand()%bestOfHowMuch;
    unsigned int winnerIndex = maleIndices[candidatIndex];
    return winnerIndex;
  }
}
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingFixedThresholdConspicific::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& /*parameters*/) const
{
  const unsigned int nMales = males.size();
  const enumSpecies femaleSpecies = female.species;
  for(unsigned int i = 0; i < nMales*2; ++i) //The female searches all males 2x
  {
    const unsigned int maleIndex = std::rand()%nMales;
    const enumSpecies maleSpecies = males[maleIndex].species;
    if (rnd::uniform() >= female.assessingError) //The female does NOT make an error
    { //The female chooses him when either having a trait smaller of bigger then 0.5
      if (femaleSpecies == piedFlycatcher)
      { //She fancies a male of speciesA
        if ( maleSpecies == piedFlycatcher) return maleIndex;
      }
      else
      { //She fancies a male of species B
        if ( maleSpecies == collaredFlycatcher) return maleIndex;
      }
    }
    else //The female makes an error
    { //The female chooses him whatever species he is
      assert(female.assessingError>0.0);
      return maleIndex;
    }
  }
  //She has not found a mate
  //Dot("Female #" + IntToStr(female.index) + " has not found a mate");
  return nMales;

}
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingFixedThresholdTraitSign::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& /*parameters*/) const
{
  const unsigned int nMales = males.size();
  const double femalePreference = female.preference;
  for(unsigned int i = 0; i < nMales*2; ++i) //The female searches all males 2x
  {
    const unsigned int maleIndex = std::rand()%nMales;
    const double maleTrait = males[maleIndex].trait;
    if (rnd::uniform() >= female.assessingError) //The female does NOT make an error
    { //The female chooses him when either having a trait smaller of bigger then 0.5
      if (femalePreference<=0.0)
      { //She fancies a male with low trait
        if ( maleTrait <= 0.0) return maleIndex;
      }
      else
      { //She fancies a male with high trait
        if ( maleTrait >= 0.0) return maleIndex;
      }
    }
    else //The female makes an error
    { //The female chooses him whatever species he is
      assert(female.assessingError>0.0);
      return maleIndex;
    }
  }
  //She has not found a mate
  return nMales;
}
//---------------------------------------------------------------------------
unsigned int StateFemaleSamplingFixedThresholdProbabilistic::getWinnerIndex(
  const std::vector<Male>& males,
  const Female& female,
  const Parameters& parameters) const
{
  const unsigned int nMales = males.size();
  const double sigmaSquared = parameters.sigmaSquared;
  const double femalePreference = female.preference;
  for(unsigned int i = 0; i < nMales*2; ++i) //The female searches all males 2x
  {
    const unsigned int maleIndex = std::rand()%nMales;
    const double maleTrait = males[maleIndex].trait;
    if (rnd::uniform() >= female.assessingError) //The female does NOT make an error
    { //The female chooses him when either having a trait smaller of bigger then 0.5
      //const double chanceToMate = std::exp(-(maleTrait-femalePreference)*(maleTrait-femalePreference)/sigmaSquared);
      const double chanceToMate = getChanceToMate(femalePreference,maleTrait,sigmaSquared);
      if (rnd::uniform() < chanceToMate) return maleIndex;
    }
    else //The female makes an error
    { //The female chooses him whatever species he is
      assert(female.assessingError>0.0);
      return maleIndex;
    }
  }
  //She has not found a mate
  return nMales;
}
//---------------------------------------------------------------------------
inline double StateFemaleSamplingFixedThresholdProbabilistic::getChanceToMate(
  const double& femalePreference,
  const double& maleTrait,
  const double& sigmaSquared)
{
   return std::exp(-(maleTrait-femalePreference)*(maleTrait-femalePreference)/sigmaSquared);
}
//---------------------------------------------------------------------------

//
// HELPER FUNCTIONS
//
//---------------------------------------------------------------------------
/*
unsigned int getBestMaleIndex(const std::vector<double>& maleTraits, const double& femalePreference, const double& assessingError)
//unsigned int getBestMaleIndex(
//  const std::vector<double>& maleTraits,
//  const double& femalePreference,
//  const double& assessingError)
{
    std::vector<unsigned int> bestMaleIndices;
    const unsigned int nCandidats = maleTraits.size();

    for (unsigned int i=0; i<nCandidats; ++i)
    {
      if (rnd::uniform() > assessingError) //The female does NOT make an error
      { //The female chooses her preference
        if ( maleTraits[i] == femalePreference) bestMaleIndices.push_back(i);
      }
      else //The female makes an error
      { //The female chooses him whatever species he is
        bestMaleIndices.push_back(i);
      }
    }
    //Select from the indices of best males the winner
    const unsigned int nCompetitors = bestMaleIndices.size();
    if (nCompetitors==0)
    {
      //Only losers, female chooses at random from all candidats
      const unsigned int bestMaleIndex = std::rand()%nCandidats;
      assert(bestMaleIndex<nCandidats);
      return bestMaleIndex;
    }
    //There are winner males, as nCompetitors > 0
    const unsigned int luckyCompetitorIndex = std::rand()%nCompetitors;
    const unsigned int bestMaleIndex = bestMaleIndices[luckyCompetitorIndex];
    assert(bestMaleIndex<nCandidats);
    return bestMaleIndex;
}
//---------------------------------------------------------------------------
//This function produces a vector of traits, in this case the select males' trait/ornaments
std::vector<double> createTraitsVector(const std::vector<Male>& males, const std::vector<unsigned int>& maleIndex)
//std::vector<double> createTraitsVector(
//  const std::vector<Male>& males,
//  const std::vector<unsigned int>& maleIndex)
{
  const unsigned int nMales = maleIndex.size();
  std::vector<double> maleTraits(nMales);
  for (unsigned int i=0; i<nMales; ++i)
  {
    const unsigned int indexMale = maleIndex[i];
    maleTraits[i] = males[indexMale].trait;
  }

  return maleTraits;
}
*/
//---------------------------------------------------------------------------
//Creates a vector of size 'bestOfHowMuch' with indices from 0 to 'nMales'
//An index might occur twice: if every index should be unique use 'createMaleIndices' instead.
std::vector<unsigned int> createIndices(
  const unsigned int& bestOfHowMuch,
  const unsigned int& nMales)
{
  std::vector<unsigned int> maleIndex(nMales);
  for(unsigned int mate = 0; mate < bestOfHowMuch; ++mate)
  {
    maleIndex[mate] = std::rand()%nMales;
    //No checking if a male is already chosen twice
    //if you want to do so, use 'createMaleIndicesUnique' instead of this function
  }
  return maleIndex;
}
//---------------------------------------------------------------------------
//Creates a vector of size 'bestOfHowMuch' from 0 to 'nMales'
//Every index is unique: if this is not important use 'createMaleIndices' instead
std::vector<unsigned int> createIndicesUnique(
  const unsigned int& bestOfHowMuch,
  const unsigned int& nMales)
{
  std::vector<unsigned int> maleIndex;

  if (bestOfHowMuch<nMales)
  {
    maleIndex.resize(bestOfHowMuch);
    assert(bestOfHowMuch<=nMales);

    maleIndex[0] = std::rand()%nMales; //0th element is always unique

    for(unsigned int mate = 1; mate < bestOfHowMuch; ++mate)
    {
      unsigned int index = 1000000; //1M
      do
      {
        index = std::rand()%nMales; //Do draw a random index
      }
      while(
        std::find(
          maleIndex.begin(),
          maleIndex.begin()+mate,
          index) != (maleIndex.begin()+mate)); //While not unqiue
      maleIndex[mate] = index;
    }
  }
  else
  { //To few males (nMales < bestOfHowMuch)
    maleIndex.resize(nMales);
    for (unsigned int male = 0; male< nMales; ++male) maleIndex[male] = male;
  }

  return maleIndex;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

helperfunctions.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef HELPERFUNCTIONS_H
#define HELPERFUNCTIONS_H
//---------------------------------------------------------------------------
#include <iostream>
#include <string>
#include <vector>
#include <sstream>
#include <assert.h>
#include <fstream>


unsigned int totalSum(const std::vector<unsigned int>& myVector);
template<class T> T totalSum(const std::vector<T>& myVector);
unsigned int getIndexLowestValue(const std::vector<double>& myVector);
unsigned int getIndexHighestValue(const std::vector<double>& myVector);
template <class T> void meanAndStdDev(const std::vector<T>& myVector, double& mean, double& stdDev);


std::string itoa(const int& x);
std::string ftoa(const double& x);
bool fileExists(const std::string& fileName);
std::string doubleToBitString(const double& anyDouble);

void Dot(const std::string& message);
void Dot(const std::vector<unsigned int>& myVector);
void Dot(const std::vector<int>& myVector);
void Dot(const std::vector<double>& myVector);



#endif // HELPERFUNCTIONS_H

 

 

 

 

 

helperfunctions.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <cmath>
#include "helperfunctions.h"
//---------------------------------------------------------------------------
//#define EXTENSIVE_LOG
//---------------------------------------------------------------------------
bool fileExists(const std::string& fileName)
{
  std::fstream fin;
  fin.open(fileName.c_str(),std::ios::in);
  if( fin.is_open() )
  {
    fin.close();
    return true;
  }
  fin.close();
  return false;
}
//---------------------------------------------------------------------------
unsigned int totalSum(const std::vector<unsigned int>& myVector)
{
  unsigned int sum=0;
  const unsigned int size = myVector.size();
  for (unsigned int i=0; i<size; ++i)
  {
    sum+=myVector[i];
  }
  return(sum);
}
//---------------------------------------------------------------------------
template<class T> T totalSum(const std::vector<T>& myVector)
{
  T sum=0;
  const unsigned int size = myVector.size();
  for (unsigned int i=0; i<size; ++i)
  {
    sum+=myVector[i];
  }
  return(sum);
}
//---------------------------------------------------------------------------
unsigned int getIndexLowestValue(const std::vector<double>& myVector)
{
  const unsigned int size = myVector.size();
  unsigned int indexLowest = 0;
  double lowestValue = myVector[0];
  for (unsigned int i=1; i<size; ++i)
  {
    const double value = myVector[i];
    if (value < lowestValue)
    {
      lowestValue = value;
      indexLowest = i;
    }
  }
  return indexLowest;
}
//---------------------------------------------------------------------------
unsigned int getIndexHighestValue(const std::vector<double>& myVector)
{
  const unsigned int size = myVector.size();
  unsigned int indexHighest = 0;
  double highestValue = myVector[0];
  for (unsigned int i=1; i<size; ++i)
  {
    const double value = myVector[i];
    if (value > highestValue)
    {
      highestValue = value;
      indexHighest = i;
    }
  }
  return indexHighest;
}
//---------------------------------------------------------------------------
std::string doubleToBitString(const double& anyDouble)
{
  static union MyUnion
  {
    unsigned int myInts[2];
    double myDouble;
  } myUnion;
  myUnion.myDouble = anyDouble;

  unsigned int myInt0 = myUnion.myInts[0];
  unsigned int myInt1 = myUnion.myInts[1];
  std::string myBitString;
  for (unsigned int i = 0; i<32; ++i)
  {
    if (myInt0%2==0)
      myBitString = "0" + myBitString;
    else
      myBitString = "1" + myBitString;
    myInt0>>=1;
  }
  for (unsigned int i = 0; i<32; ++i)
  {
    if (myInt1%2==0)
      myBitString = "0" + myBitString;
    else
      myBitString = "1" + myBitString;
    myInt1>>=1;
  }
  return myBitString;
}
//---------------------------------------------------------------------------
template <class T>
void meanAndStdDev(const std::vector<T>& myVector, double& mean, double& stdDev)
{
  const unsigned int size = myVector.size();
  assert(size>1);
  const double dSize = static_cast<double>(size);
  mean = 0;
  double sumX = 0;
  double sumXsquared = 0;
  for (unsigned int i=0; i<size; ++i)
  {
    const T value = myVector[i];
    sumX+=value;
    sumXsquared+=(value*value);
    mean+=value;
  }

  mean/=dSize;
  stdDev = std::sqrt(((dSize*sumXsquared)-(sumX*sumX))/(dSize *(dSize-1.0)));
}
//---------------------------------------------------------------------------
std::string ftoa(const double& x)
{
  std::ostringstream o;
  if (!(o << x)) return "ERROR";
  return o.str();
}
//---------------------------------------------------------------------------
std::string itoa(const int& x)
{
  std::ostringstream o;
  if (!(o << x)) return "ERROR";
  return o.str();
}
//---------------------------------------------------------------------------
void Dot(const std::string& message)
{
  #ifndef NDEBUG
    std::clog << "ODS: " << message << std::endl;
  #endif
}
//---------------------------------------------------------------------------
void Dot(const std::vector<double>& myVector)
{
  #ifndef NDEBUG
  std::string output = "index: ";
  const unsigned int size = myVector.size();
  for (unsigned int i=0; i<size; ++i)
  {
    output+=ftoa(myVector[i])+" ";
  }
  Dot(output);
  #endif
}
//---------------------------------------------------------------------------
void Dot(const std::vector<int>& myVector)
{
  #ifndef NDEBUG
  std::string output = "index: ";
  const unsigned int size = myVector.size();
  for (unsigned int i=0; i<size; ++i)
  {
    output+=itoa(myVector[i])+" ";
  }
  Dot(output);
  #endif
}
//---------------------------------------------------------------------------
void Dot(const std::vector<unsigned int>& myVector)
{
  #ifndef NDEBUG
  std::string output = "index: ";
  const unsigned int size = myVector.size();
  for (unsigned int i=0; i<size; ++i)
  {
    output+=itoa(myVector[i])+" ";
  }
  Dot(output);
  #endif
}
//---------------------------------------------------------------------------

 

 

 

 

 

main.cpp

 

#include <iostream>
#include <memory>
#include <vector>
#include <fstream>
#include <string>
#include <assert.h>

#include "parameters.h"
#include "simulation.h"
#include "helperfunctions.h"
//Parameters readParametersFromFile(const std::string& fileName);
//---------------------------------------------------------------------------
int main(int argc, char* argv[])
{
  if (argc<2)
  {
    std::cout << "Please enter name of data file after file name,"
      << "\ne.g. ProjectChrisWileySTL data.txt" << std::endl;
      return 0;
  }
  if (fileExists(argv[1])==false)
  {
    std::cout << "Please enter a VALID and EXISTING name of data file after file name." << std::endl;
    return 0;
  }

  //Parameters parameters = readParametersFromFile(argv[1]);

  //Simulation simulation(parameters);
  //std::auto_ptr<Simulation> simulation(new Simulation(
}
//---------------------------------------------------------------------------

 

 

 

 

 

matetally.h

 

#ifndef MATETALLY_H
#define MATETALLY_H

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <string>
#include <vector>
#include "bird.h"
//---------------------------------------------------------------------------
class MateTally
{
  public:
  MateTally() { reset(); }
  void tally(const Bird& mother, const Bird& father);
  void reset();
  std::string get() const;
  bool isNull() const;
  double calculateBiasA() const;
  double calculateBiasB() const;
  double calculateFractionMixedPairs() const;

  unsigned int getNmateAA() const { return nMateAA; } //Female A - male A
  unsigned int getNmateAB() const { return nMateAB; } //Female A - male B
  unsigned int getNmateBA() const { return nMateBA; } //Female B - male A
  unsigned int getNmateBB() const { return nMateBB; } //Female B - male B
  unsigned int getNmateAll() const { return nMateAA + nMateAB + nMateBA + nMateBB; }

  void operator/=(const unsigned int& intValue);
  void operator+=(const MateTally&  mateTally);

  private:
  unsigned int nMateAA; //Female A - male A
  unsigned int nMateAB; //Female A - male B
  unsigned int nMateBA; //Female B - male A
  unsigned int nMateBB; //Female B - male B
};
//---------------------------------------------------------------------------
#endif // MATETALLY_H

 

 

 

 

 

matetally.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "matetally.h"
//---------------------------------------------------------------------------
void MateTally::reset()
{
  nMateAA = 0;
  nMateAB = 0;
  nMateBA = 0;
  nMateBB = 0;
}
//---------------------------------------------------------------------------
void MateTally::operator/=(const unsigned int& intValue)
{
  nMateAA/=intValue;
  nMateAB/=intValue;
  nMateBA/=intValue;
  nMateBB/=intValue;
}
//---------------------------------------------------------------------------
void MateTally::operator+=(const MateTally& mateTally)
{
  nMateAA+=mateTally.nMateAA;
  nMateAB+=mateTally.nMateAB;
  nMateBA+=mateTally.nMateBA;
  nMateBB+=mateTally.nMateBB;
}
//---------------------------------------------------------------------------
std::string MateTally::get() const
{
  std::string result;
  result+= "nMateAA: " + itoa(nMateAA);
  result+= " ,nMateAB: " + itoa(nMateAB);
  result+= " ,nMateBA: " + itoa(nMateBA);
  result+= " ,nMateBB: " + itoa(nMateBB);
  return result;
}
//---------------------------------------------------------------------------
bool MateTally::isNull() const
{
  if ( nMateAA!=0
    || nMateAB!=0
    || nMateBA!=0
    || nMateBB!=0) return false;
  return true;
}
//---------------------------------------------------------------------------
void MateTally::tally(const Bird& mother, const Bird& father)
{
  if (mother.species == piedFlycatcher)
  {
    if (father.species == piedFlycatcher)
      { ++nMateAA; }
    else
      { ++nMateAB; }
  }
  else
  {
    if (father.species == piedFlycatcher)
    { ++nMateBA; }
    else
    { ++nMateBB; }
  }
}
//---------------------------------------------------------------------------
double MateTally::calculateBiasA() const
{
  return (nMateAB + nMateBA == 0 ? 0.0
    : static_cast<double>(nMateAB) / static_cast<double>(nMateAB + nMateBA));

}
//---------------------------------------------------------------------------
double MateTally::calculateBiasB() const
{
  return (nMateAB + nMateBA == 0 ? 0.0
    : static_cast<double>(nMateBA) / static_cast<double>(nMateAB + nMateBA));
}
//---------------------------------------------------------------------------
double MateTally::calculateFractionMixedPairs() const
{
  return (this->isNull()==true ? 0.0
    : static_cast<double>(nMateAB + nMateBA) / static_cast<double>(nMateAA + nMateAB + nMateBA + nMateBB));
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

matingsystem.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef MATINGSYSTEM_H
#define MATINGSYSTEM_H
//---------------------------------------------------------------------------
#include <iostream>
#include <vector>
#include <memory>
#include <string>
#include <algorithm>


#include "bird.h"
#include "parameters.h"
#include "femalesampling.h"
#include "matetally.h"
#include "helperfunctions.h"
#include "random.h"

class StateMatingSystemBase
{
  public:
  StateMatingSystemBase() {}
  virtual std::string getMatingSystem() const = 0;
  virtual Offspring mate(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
    MateTally& tally
    ) = 0;
};
//---------------------------------------------------------------------------
class StateMatingSystemMonogamyFixedNumberOffspring : public StateMatingSystemBase
{
  public:
  StateMatingSystemMonogamyFixedNumberOffspring() {}
  std::string getMatingSystem() const { return "MonogamyFixedNumberOffspring"; }
  Offspring mate(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
    MateTally& tally
    );

};
//---------------------------------------------------------------------------
class StateMatingSystemPolygynyFixedNumberOffspring : public StateMatingSystemBase
{
  public:
  StateMatingSystemPolygynyFixedNumberOffspring() {}
  std::string getMatingSystem() const { return "PolygynyFixedNumberOffspring"; }
  Offspring mate(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
    MateTally& tally
    );
};
//---------------------------------------------------------------------------
class StateMatingSystemMonogamyFreeNumberOffspring : public StateMatingSystemBase
{
  public:
  StateMatingSystemMonogamyFreeNumberOffspring() {}
  std::string getMatingSystem() const { return "MonogamyFreeNumberOffspring"; }
  Offspring mate(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
    MateTally& tally
    );

};
//---------------------------------------------------------------------------
class StateMatingSystemPolygynyFreeNumberOffspring : public StateMatingSystemBase
{
  public:
  StateMatingSystemPolygynyFreeNumberOffspring() {}
  std::string getMatingSystem() const { return "PolygynyFreeNumberOffspring"; }
  Offspring mate(
    std::vector<Female>& females,
    std::vector<Male>& males,
    const Parameters& parameters,
    const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
    MateTally& tally
    );
};
//---------------------------------------------------------------------------
#endif // MATINGSYSTEM_H

 

 

 

 

 

matingsystem.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "matingsystem.h"
#include "speciestally.h"
//---------------------------------------------------------------------------
Offspring StateMatingSystemMonogamyFixedNumberOffspring::mate(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters,
  const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
  MateTally& mateTally
  )
{
  std::auto_ptr<MateTally> debugMateTally(new MateTally);

  assert(mateTally.isNull()==true);
  //Get the fraction of species A and B
  const SpeciesTally speciesTally(males,females);

  unsigned int nMalesAwanted   = parameters.getNmalesAwanted();
  unsigned int nMalesBwanted   = parameters.getNmalesBwanted();
  assert(nMalesAwanted + nMalesBwanted == parameters.nMales);
  unsigned int nFemalesAwanted = parameters.getNfemalesAwanted();
  unsigned int nFemalesBwanted = parameters.getNfemalesBwanted();
  assert(nFemalesAwanted + nFemalesBwanted == parameters.nFemales);

  const unsigned int nSpeciesAwanted = nMalesAwanted + nFemalesAwanted;
  const unsigned int nSpeciesBwanted = nMalesBwanted + nFemalesBwanted;
  const double mutation = parameters.mutationRate;

  assert(females.size()>0);

  //Create the couples
  std::vector< Couple > couples;
  for(unsigned int female = 0; female < females.size() && males.size() > 0; /*nothing*/)
  {
    //Let the female pick her favorite
    const unsigned int winnerIndex = femaleSampling->getWinnerIndex(males, females[female], parameters);
    assert(winnerIndex <= males.size());
    //Dot("The index of the mate in vector 'males' is: " + IntToStr(winnerIndex));
    if (winnerIndex == males.size())
    {
      //The female did not mate
      assert(parameters.femaleSampling!=bestOfNconspicific && parameters.femaleSampling!=bestOfNextremeTrait && parameters.femaleSampling!=bestOfNclosestTrait);
      ++female;
      continue;
    }
    //Form the std::pair from the lovers and store them in the std::vector couples
    Couple couple(females[female],males[winnerIndex]);
    couples.push_back(couple);
    //Remove the couple from the desperate others
    males.erase(males.begin() + winnerIndex);
    females.erase(females.begin() + female);
  }

  //Let the couples produce offspring
  std::vector<Bird> speciesA; //speciesA.reserve(nSpeciesAwanted);
  std::vector<Bird> speciesB; //speciesB.reserve(nSpeciesBwanted);
  const unsigned int nCouples = couples.size();
  const unsigned int maxTries = (nSpeciesAwanted + nSpeciesBwanted) * (nSpeciesAwanted + nSpeciesBwanted);
  unsigned int myTry = 0;
  while(speciesA.size()!=nSpeciesAwanted || speciesB.size()!=nSpeciesBwanted)
  {
    ++myTry;
    if (myTry==maxTries) break;
    for(unsigned int couple = 0; couple < nCouples; ++couple)
    {
      //Find out whether he is a conspicific and put it in the results
      Bird birdy(couples[couple].female,couples[couple].male,mutation);
      if (birdy.species == piedFlycatcher)
      { //Birdy is of speciesA
        debugMateTally->tally(couples[couple].female,couples[couple].male);
        if (speciesA.size() < nSpeciesAwanted)
        { //When we can use offspring of speciesA
          speciesA.push_back(birdy);
          mateTally.tally(couples[couple].female,couples[couple].male);
        }
      }
      else
      { //Birdy is of speciesB
        debugMateTally->tally(couples[couple].female,couples[couple].male); //DEBUG
        if (speciesB.size() < nSpeciesBwanted)
        { //When we can use offspring of speciesA
          speciesB.push_back(birdy);
          mateTally.tally(couples[couple].female,couples[couple].male);
        }
      }
    }
  }

  Dot("Broke while loop in reproduction");
  //Now the vectors of speciesA and speciesB are created,
  //they have to be changed to vectors of males and females
  Offspring offspring;
  if (myTry==maxTries)
  {
    //Create offspring of a failed simulation.
    //Easy! Just give the empty Offspring back
    return offspring;
  }

  assert(speciesA.size()==nSpeciesAwanted);
  assert(speciesB.size()==nSpeciesBwanted);
  //Females
  for (unsigned int i=0; i<nFemalesAwanted; ++i) offspring.females.push_back(speciesA[i]);
  for (unsigned int i=0; i<nFemalesBwanted; ++i) offspring.females.push_back(speciesB[i]);
  //Males
  for (unsigned int i=nFemalesAwanted; i<nSpeciesAwanted; ++i) offspring.males.push_back(speciesA[i]);
  for (unsigned int i=nFemalesBwanted; i<nSpeciesBwanted; ++i) offspring.males.push_back(speciesB[i]);

  assert(offspring.females.size()==nFemalesAwanted+nFemalesBwanted);
  assert(offspring.males.size()==nMalesAwanted+nMalesBwanted);

  //Dot("mateTally: "+mateTally.get());
  //Dot("DebugTally: "+debugMateTally->get());
  return offspring;

}
//---------------------------------------------------------------------------
Offspring StateMatingSystemPolygynyFixedNumberOffspring::mate(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters,
  const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
  MateTally& mateTally
  )
{
  assert(mateTally.isNull()==true);
  //Get the fraction of species A and B
  const SpeciesTally speciesTally(males,females);
  /*
  const unsigned int nMalesAwanted   = static_cast<double>(parameters.nMales * speciesTally.getNallMalesA()) / static_cast<double>(speciesTally.getNallMales());
  const unsigned int nMalesBwanted   = static_cast<double>(parameters.nMales * speciesTally.getNallMalesB()) / static_cast<double>(speciesTally.getNallMales());
  const unsigned int nFemalesAwanted = static_cast<double>(parameters.nFemales * speciesTally.getNfemalesA()) / static_cast<double>(speciesTally.getNallFemales());
  const unsigned int nFemalesBwanted = static_cast<double>(parameters.nFemales * speciesTally.getNfemalesB()) / static_cast<double>(speciesTally.getNallFemales());
  */
  unsigned int nMalesAwanted   = parameters.getNmalesAwanted();
  unsigned int nMalesBwanted   = parameters.getNmalesBwanted();
  assert(nMalesAwanted + nMalesBwanted == parameters.nMales);
  unsigned int nFemalesAwanted = parameters.getNfemalesAwanted();
  unsigned int nFemalesBwanted = parameters.getNfemalesBwanted();
  assert(nFemalesAwanted + nFemalesBwanted == parameters.nFemales);

  const unsigned int nSpeciesAwanted = nMalesAwanted + nFemalesAwanted;
  const unsigned int nSpeciesBwanted = nMalesBwanted + nFemalesBwanted;
  const double mutation = parameters.mutationRate;
  const unsigned int nFemales = females.size();
  assert(nFemales>0);
  std::vector<Bird> speciesA;
  std::vector<Bird> speciesB;

  const unsigned int maxTries = (nSpeciesAwanted + nSpeciesBwanted) * (nSpeciesAwanted + nSpeciesBwanted);
  unsigned int myTry = 0;
  while(speciesA.size()!=nSpeciesAwanted || speciesB.size()!=nSpeciesBwanted)
  {
    ++myTry;
    if (myTry==maxTries) break;
    for(unsigned int female = 0; female < nFemales; ++female)
    {
      //Let the female pick her favorite
      const unsigned int winnerIndex = femaleSampling->getWinnerIndex(males, females[female], parameters);
      assert(winnerIndex <= males.size());

      //Dot("The index of the mate in vector 'males' is: " + IntToStr(winnerIndex));
      if (winnerIndex == males.size())
      {
        //The female did not mate
        assert(parameters.femaleSampling!=bestOfNconspicific && parameters.femaleSampling!=bestOfNextremeTrait && parameters.femaleSampling!=bestOfNclosestTrait);
        continue;
      }

      //Find out whether he is a conspicific and put it in the results
      Bird birdy(females[female],males[winnerIndex],mutation);
      if (birdy.species == piedFlycatcher)
      { //Birdy is of speciesA
        if (speciesA.size() < nSpeciesAwanted)
        { //When we can use offspring of speciesA
          speciesA.push_back(birdy);
          mateTally.tally(females[female],males[winnerIndex]);
        }
      }
      else
      { //Birdy is of speciesB
        if (speciesB.size() < nSpeciesBwanted)
        { //When we can use offspring of speciesA
          speciesB.push_back(birdy);
          mateTally.tally(females[female],males[winnerIndex]);
        }
      }
    }
  }

  Dot("Broke while loop in reproduction above line ");

  //Now the vectors of speciesA and speciesB are created,
  //they have to be changed to vectors of males and females

  Offspring offspring;
  if (myTry==maxTries)
  {
    //Create offspring of a failed simulation.
    //Easy! Just give the empty Offspring back
    return offspring;
  }

  assert(speciesA.size()==nSpeciesAwanted);
  assert(speciesB.size()==nSpeciesBwanted);
  //Females
  for (unsigned int i=0; i<nFemalesAwanted; ++i) offspring.females.push_back(speciesA[i]);
  for (unsigned int i=0; i<nFemalesBwanted; ++i) offspring.females.push_back(speciesB[i]);
  //Males
  for (unsigned int i=nFemalesAwanted; i<nSpeciesAwanted; ++i) offspring.males.push_back(speciesA[i]);
  for (unsigned int i=nFemalesBwanted; i<nSpeciesBwanted; ++i) offspring.males.push_back(speciesB[i]);

  assert(offspring.females.size()==nFemalesAwanted+nFemalesBwanted);
  assert(offspring.males.size()==nMalesAwanted+nMalesBwanted);
  return offspring;
}
//---------------------------------------------------------------------------
Offspring StateMatingSystemMonogamyFreeNumberOffspring::mate(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters,
  const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
  MateTally& mateTally
  )
{
  std::auto_ptr<MateTally> debugMateTally(new MateTally);

  assert(mateTally.isNull()==true);
  const double mutation = parameters.mutationRate;
  assert(females.size()>0);

  //Create the couples
  std::vector< Couple > couples;
  for(unsigned int female = 0; female < females.size() && males.size() > 0; )
  {
    //Let the female pick her favorite
    const unsigned int winnerIndex = femaleSampling->getWinnerIndex(males, females[female], parameters);
    assert(winnerIndex <= males.size());
    //Dot("The index of the mate in vector 'males' is: " + IntToStr(winnerIndex));
    if (winnerIndex == males.size())
    {
      //The female did not mate
      assert(parameters.femaleSampling!=bestOfNconspicific && parameters.femaleSampling!=bestOfNextremeTrait && parameters.femaleSampling!=bestOfNclosestTrait);
      ++female;
      continue;
    }
    //Form the std::pair from the lovers and store them in the std::vector couples
    Couple couple(females[female],males[winnerIndex]);
    couples.push_back(couple);
    //Remove the couple from the desperate others
    males.erase(males.begin() + winnerIndex);
    females.erase(females.begin() + female);
  }

  //Let the couples produce offspring
  Offspring offspring;
  const unsigned int nCouples = couples.size();
  for(unsigned int couple = 0; couple < nCouples; ++couple)
  {
    const unsigned int nOffspring = parameters.nOffspring;
    for (unsigned i=0; i<nOffspring; ++i)
    {
      //Find out whether he is a conspicific and put it in the results
      Bird birdy(couples[couple].female,couples[couple].male,mutation);
      debugMateTally->tally(couples[couple].female,couples[couple].male); //DEBUG
      mateTally.tally(couples[couple].female,couples[couple].male);
      if (std::rand()%2==0)
        offspring.males.push_back(birdy);
      else
        offspring.females.push_back(birdy);
    }
  }

  //Now the vectors of speciesA and speciesB are created,
  //they have to be changed to vectors of males and females
  //Dot("mateTally: "+mateTally.get());
  //Dot("DebugTally: "+debugMateTally->get());
  return offspring;
}
//---------------------------------------------------------------------------
Offspring StateMatingSystemPolygynyFreeNumberOffspring::mate(
  std::vector<Female>& females,
  std::vector<Male>& males,
  const Parameters& parameters,
  const std::auto_ptr<StateFemaleSamplingBase>& femaleSampling,
  MateTally& mateTally
  )
{
  std::auto_ptr<MateTally> debugMateTally(new MateTally);
  assert(mateTally.isNull()==true);
  //Get the fraction of species A and B
  const double mutation = parameters.mutationRate;
  const unsigned int nFemales = females.size();
  assert(nFemales>0);
  Offspring offspring;
  for(unsigned int female = 0; female < nFemales; ++female)
  {
    //Let the female pick her favorite
    const unsigned int winnerIndex = femaleSampling->getWinnerIndex(males, females[female], parameters);
    assert(winnerIndex <= males.size());

    //Dot("The index of the mate in vector 'males' is: " + IntToStr(winnerIndex));
    if (winnerIndex == males.size())
    {
      //The female did not mate
      assert(parameters.femaleSampling!=bestOfNconspicific && parameters.femaleSampling!=bestOfNextremeTrait && parameters.femaleSampling!=bestOfNclosestTrait);
      continue;
    }

    const unsigned int nOffspring = parameters.nOffspring;
    //Let them produce many offspring
    for (unsigned i=0; i<nOffspring; ++i)
    {
      //Find out whether he is a conspicific and put it in the results
      Bird birdy(females[female],males[winnerIndex],mutation);
      debugMateTally->tally(females[female],males[winnerIndex]); //DEBUG
      mateTally.tally(females[female],males[winnerIndex]);
      if (std::rand()%2==0)
        offspring.males.push_back(birdy);
      else
        offspring.females.push_back(birdy);
    }
  }

  return offspring;
}
//---------------------------------------------------------------------------

 

 

 

 

 

parameters.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef PARAMETERS_H
#define PARAMETERS_H
//---------------------------------------------------------------------------
#include "enums.h"
#include "helperfunctions.h"

class Parameters
{
  public:
  unsigned int nGenerations;  //Number of generations
  unsigned int nFemales;      //Number of females
  unsigned int nMales;        //Number of males
  double fractionMaleA;       //Fraction of males being of species A
  double fractionFemaleA;     //Fraction of males being of species B
  unsigned int bestOfHowMuch; //Value of N in Best-of-N
  double assessingErrorA ;    //The probability a female A assess a male wrong
  double assessingErrorB;     //The probability a female B assess a male wrong
  double meanTraitA;
  double meanTraitB;
  double meanPreferenceA;
  double meanPreferenceB;
  double stdDevTraitA;
  double stdDevTraitB;
  double stdDevPreferenceA;
  double stdDevPreferenceB;
  double sigmaSquared;
  unsigned int nSimulations; //Number of simulations
  unsigned int nOffspring;   //Number of offspring produced per couple
  //unsigned int simIndex; //The simulation's index, or: the simIndex-th simulation that is being run
  double surviveSpeciesAlpha;
  double surviveSpeciesBeta;
  double costTrait;
  double costPreference;
  double mutationRate;
  enumMatingSystem matingSystem;     //Monogamy/polygyny
  enumFemaleSampling femaleSampling; //The way a female samples the male population
  enumDensityDependentSelection densityDependentSelection;
  void reset();
  unsigned int getNmalesAwanted() const;
  unsigned int getNmalesBwanted() const;
  unsigned int getNfemalesAwanted() const;
  unsigned int getNfemalesBwanted() const;
  void readFromFile(const std::string& fileName);
  void writeToFile(const std::string& fileName);
};




#endif // PARAMETERS_H

 

 

 

 

 

parameters.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include <cstdlib>
//---------------------------------------------------------------------------
#include "parameters.h"
//---------------------------------------------------------------------------
void Parameters::reset()
{
   nFemales = 0;
   nMales = 0;
   fractionMaleA = 0.0;
   fractionFemaleA = 0.0;
   bestOfHowMuch = 0;
   assessingErrorA = 0.0;
   assessingErrorB = 0.0;
   meanTraitA = 0.0;
   meanTraitB = 0.0;
   meanPreferenceA = 0.0;
   meanPreferenceB = 0.0;
   stdDevTraitA = 0.0;
   stdDevTraitB = 0.0;
   stdDevPreferenceA = 0.0;
   stdDevPreferenceB = 0.0;
   sigmaSquared = 0.0;
   nSimulations = 0;
   //simIndex = 0;
   costTrait = 0.0;
   costPreference = 0.0;
   mutationRate = 0.0;
   matingSystem = monogamy;
   femaleSampling = bestOfNconspicific;
   densityDependentSelection = afterMating;
}
//---------------------------------------------------------------------------
unsigned int Parameters::getNmalesAwanted() const
{
  //return static_cast<double>(nMales) * fractionMaleA;
  //return nMales * fractionMaleA;
  const unsigned int result = nMales * fractionMaleA;
  const unsigned int nMalesBwanted = getNmalesBwanted();
  if (result+nMalesBwanted==nMales) return result;
  else if (result+nMalesBwanted==nMales-1) return result+1;
  else assert(!"Should not get here"); std::exit(1); return 0;
}
//---------------------------------------------------------------------------
unsigned int Parameters::getNmalesBwanted() const
{
  return nMales * (1.0-fractionMaleA);
}
//---------------------------------------------------------------------------
unsigned int Parameters::getNfemalesAwanted() const
{
  //return static_cast<double>(nFemales) * fractionFemaleA;
  //return nFemales * fractionFemaleA;
  const unsigned int result = nFemales * fractionFemaleA;
  const unsigned int nFemalesBwanted = getNfemalesBwanted();
  if (result+nFemalesBwanted==nFemales) return result;
  else if (result+nFemalesBwanted==nFemales-1) return result+1;
  else assert(!"Should not get here"); std::exit(1); return 0;
}
//---------------------------------------------------------------------------
unsigned int Parameters::getNfemalesBwanted() const
{
  //return static_cast<double>(nFemales) * (1.0-fractionFemaleA);
  return nFemales * (1.0-fractionFemaleA);
}
//---------------------------------------------------------------------------
void Parameters::readFromFile(const std::string& fileName)
{
  assert(fileExists(fileName)==true);
  std::ifstream in (fileName.c_str());
  std::string myString;
  int tempInt;
  for (int i=0; !in.eof(); ++i)
  {
    in >> myString;
    if (myString=="nGenerations")      { in >> nGenerations; continue; }
    if (myString=="nFemales")          { in >> nFemales; continue; }
    if (myString=="nMales")            { in >> nMales; continue; }
    if (myString=="fractionMaleA")     { in >> fractionMaleA; continue; }
    if (myString=="fractionFemaleA")   { in >> fractionFemaleA; continue; }
    if (myString=="bestOfHowMuch")     { in >> bestOfHowMuch; continue; }
    if (myString=="assessingErrorA")   { in >> assessingErrorA; continue; }
    if (myString=="assessingErrorB")   { in >> assessingErrorB; continue; }
    if (myString=="meanTraitA")        { in >> meanTraitA; continue; }
    if (myString=="meanTraitB")        { in >> meanTraitB; continue; }
    if (myString=="meanPreferenceA")   { in >> meanPreferenceA; continue; }
    if (myString=="meanPreferenceB")   { in >> meanPreferenceB; continue; }
    if (myString=="stdDevTraitA")      { in >> stdDevTraitA; continue; }
    if (myString=="stdDevTraitB")      { in >> stdDevTraitB; continue; }
    if (myString=="stdDevPreferenceA") { in >> stdDevPreferenceA; continue; }
    if (myString=="stdDevPreferenceB") { in >> stdDevPreferenceB; continue; }
    if (myString=="sigmaSquared")      { in >> sigmaSquared; continue; }
    if (myString=="nSimulations")      { in >> nSimulations; continue; }
    if (myString=="nOffspring")        { in >> nOffspring; continue; }
    if (myString=="surviveSpeciesAlpha") { in >> surviveSpeciesAlpha; continue; }
    if (myString=="surviveSpeciesBeta")  { in >> surviveSpeciesBeta; continue; }
    if (myString=="costTrait")  { in >> costTrait; continue; }
    if (myString=="costPreference")  { in >> costPreference; continue; }
    if (myString=="mutationRate")  { in >> mutationRate; continue; }
    if (myString=="matingSystem")  { in >> tempInt; matingSystem = static_cast<enumMatingSystem>(tempInt); continue; }
    if (myString=="femaleSampling")  { in >> tempInt; femaleSampling = static_cast<enumFemaleSampling>(tempInt); continue; }
    if (myString=="densityDependentSelection")  { in >> tempInt; densityDependentSelection = static_cast<enumDensityDependentSelection>(tempInt); continue; }
    assert(!"Unknown file parameter"); std::exit(1);
  }
  in.close();
}
//---------------------------------------------------------------------------
void Parameters::writeToFile(const std::string& fileName)
{
  std::ofstream out (fileName.c_str());
  out << "nGenerations" << '\t'       << nGenerations << '\n';
  out << "nFemales" << '\t'           << nFemales << '\n';
  out << "nMales" << '\t'             << nMales << '\n';
  out << "fractionMaleA" << '\t'      << fractionMaleA << '\n';
  out << "fractionFemaleA" << '\t'    << fractionFemaleA << '\n';
  out << "bestOfHowMuch" << '\t'      << bestOfHowMuch << '\n';
  out << "assessingErrorA" << '\t'    << assessingErrorA << '\n';
  out << "assessingErrorB" << '\t'    << assessingErrorB << '\n';
  out << "meanTraitA" << '\t'         << meanTraitA << '\n';
  out << "meanTraitB" << '\t'         << meanTraitB << '\n';
  out << "meanPreferenceA" << '\t'    << meanPreferenceA << '\n';
  out << "meanPreferenceB" << '\t'    << meanPreferenceB << '\n';
  out << "stdDevTraitA" << '\t'       << stdDevTraitA << '\n';
  out << "stdDevTraitB" << '\t'       << stdDevTraitB << '\n';
  out << "stdDevPreferenceA" << '\t'  << stdDevPreferenceA << '\n';
  out << "stdDevPreferenceB" << '\t'  << stdDevPreferenceB << '\n';
  out << "sigmaSquared" << '\t'       << sigmaSquared << '\n';
  out << "nSimulations" << '\t'       << nSimulations << '\n';
  out << "nOffspring" << '\t'         << nOffspring << '\n';
  out << "surviveSpeciesAlpha" << '\t'  << surviveSpeciesAlpha << '\n';
  out << "surviveSpeciesBeta" << '\t'   << surviveSpeciesBeta << '\n';
  out << "costTrait" << '\t'   << costTrait << '\n';
  out << "costPreference" << '\t'   << costPreference << '\n';
  out << "mutationRate" << '\t'   << mutationRate << '\n';
  out << "matingSystem" << '\t'   << static_cast<int>(matingSystem) << '\n';
  out << "femaleSampling" << '\t'   << static_cast<int>(femaleSampling) << '\n';
  out << "densityDependentSelection" << '\t'   << static_cast<int>(densityDependentSelection) << '\n';
  out.close();
}
//---------------------------------------------------------------------------

 

 

 

 

 

random.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef PROJECTCHRISWILEY_RANDOM_H
#define PROJECTCHRISWILEY_RANDOM_H
//---------------------------------------------------------------------------
#include <cmath>
#include <cstdlib>
//---------------------------------------------------------------------------
class rnd
{
  private:
  static long int idum;
  static int MBIG;
  static int MSEED;
  static int MZ;
  static double FAC;


  public:
  //Returns a uniform random value from 0.0 to 1.0
  static double uniform()
  {
    static int inext,inextp;
    static long ma[56];
    static int iff=0;
    long mj,mk;
    int i,ii,k;
    if (idum < 0 || iff == 0)
      {
      iff=1;
      mj=MSEED-(idum < 0 ? -idum : idum);
      mj %= MBIG;
      ma[55]=mj;
      mk=1;
      for (i=1;i<=54;i++)
        {
        ii=(21*i) % 55;
        ma[ii]=mk;
        mk=mj-mk;
        if (mk < MZ) mk += MBIG;
        mj=ma[ii];
        }
      for (k=1;k<=4;k++)
        for (i=1;i<=55;i++)
          {
          ma[i] -= ma[1+(i+30) % 55];
          if (ma[i] < MZ) ma[i] += MBIG;
          }
      inext=0;
      inextp=31;
      idum=1;
      }
    if (++inext == 56) inext=1;
    if (++inextp == 56) inextp=1;
    mj=ma[inext]-ma[inextp];
    if (mj < MZ) mj += MBIG;
    ma[inext]=mj;
    return mj*FAC;
  }

  static double gasdev()
  {
    static int iset=0;
    static double gset;
    double fac,r,v1,v2;
    //double uniform();
    uniform();
    if  (iset == 0)
    {
      do
      {
        v1=2.0*uniform()-1.0;
        v2=2.0*uniform()-1.0;
        r=v1*v1+v2*v2;
      }
      while (r >= 1.0);
      fac=sqrt(-2.0*log(r)/r);
      gset=v1*fac;
      iset=1;
      return v2*fac;
    }
    else
    {
      iset=0;
      return gset;
    }
  }
  //Sets the seed for the random number sequence
  static void SetSeed(const int& seed)
  {
    int i;
    idum = long(-std::abs(seed));
    for (i=0; i<100; i++) uniform();
  }
  //Returns a random number from a gaussian distrubution
  static double normal(const double& mean,const double& stdev)
  {
    return gasdev()*stdev + mean;
  }
  //Returns a random integer from 0 to N
  static int RandomNumber(const int& N)
  {
    double x;
    int out;
    x=std::floor(uniform()*N);
    out=int (x);
    return out;
  }
};


#endif // PROJECTCHRISWILEY_RANDOM_H

 

 

 

 

 

random.cpp

 

#include "random.h"

long int rnd::idum;
int rnd::MBIG  = 1000000000;
int rnd::MSEED = 161803398;
int rnd::MZ = 0;
double rnd::FAC = (1.0/MBIG);

 

 

 

 

 

results.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef RESULTS_H
#define RESULTS_H
//---------------------------------------------------------------------------
#include <vector>
//---------------------------------------------------------------------------
class ResultsSingleton
{
  public:
  static ResultsSingleton* instance()
  {
    if (mpInstance==0) mpInstance = new ResultsSingleton();
    return mpInstance;
  }


  protected:
  ResultsSingleton() {}
  private:
  static ResultsSingleton* mpInstance;
};



#endif // RESULTS_H

 

 

 

 

 

results.cpp

 

#include "results.h"

 

 

 

 

 

simulation.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef PROJECTCHRISWILEY_SIMULATION_H
#define PROJECTCHRISWILEY_SIMULATION_H
//---------------------------------------------------------------------------
#include <vector>
#include <memory>
#include <cmath>
#include <string>
#include <algorithm>
#include "enums.h"
#include "parameters.h"
#include "matetally.h"
#include "speciestally.h"
#include "helperfunctions.h"
#include "timeseries.h"
#include "bird.h"
//---------------------------------------------------------------------------
//Forward declarations of States
class StateMatingSystemBase;
class StateFemaleSamplingBase;
class StateDensityDependentSelectionBase;

class Simulation
{
  public:
  Simulation(const Parameters&);
  void execute();
  TimeSeries getTimeSeries() const { return mTimeSeries; }
  void setMatingSystem(const enumMatingSystem&, const enumDensityDependentSelection&);

  void setFemaleSampling(const enumFemaleSampling&);
  void setDensityDependentSelection(const enumDensityDependentSelection&);

  //void showPopulation(TStringGrid * stringGrid) const;
  static inline double chanceToSurviveSpecies(const double& descent, const double& alpha, const double& beta);
  static inline double chanceToSurviveTrait(const double& trait, const double& costTrait);
  static inline double chanceToSurvivePreference(const double& preference, const double& costPreference);

  //Debugging
  private:
  const Parameters mParameters;

  void speciesSelection();
  void viabilitySelection();
  bool willDieSpecies(const Bird& bird);
  bool willDiePreference(const Female& female);
  bool willDieTrait(const Male& male);

  void mating();
  std::auto_ptr<StateMatingSystemBase> mMatingSystem;
  std::auto_ptr<StateFemaleSamplingBase> mFemaleSampling;
  std::auto_ptr<StateDensityDependentSelectionBase> mDensityDependentSelection;
  std::auto_ptr<MateTally> mMateTally;
  std::vector<Bird> mMales;
  std::vector<Bird> mFemales;
  std::vector<Bird> mOffspring;
  TimeSeries mTimeSeries;

  public:
  std::string mError;
};

#include "femalesampling.h"
#include "matingsystem.h"
#include "densitydependentselection.h"

#endif // PROJECTCHRISWILEY_SIMULATION_H

 

 

 

 

 

simulation.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "bird.h"
#include "simulation.h"
//---------------------------------------------------------------------------
Simulation::Simulation(const Parameters& parameters)
  : mParameters(parameters),
  mMateTally(new MateTally),
  mError("Finished simulation without errors.")
{
  //Copy parameters from mParameters
  setMatingSystem(mParameters.matingSystem, mParameters.densityDependentSelection);
  setFemaleSampling(mParameters.femaleSampling);
  //Reset results
  mMateTally->reset();
  mMales = Bird::createMales(mParameters);
  mFemales = Bird::createFemales(mParameters);

  //Shuffle the males and females
  std::random_shuffle(mMales.begin()  , mMales.end()  );
  std::random_shuffle(mFemales.begin(), mFemales.end());
}
//---------------------------------------------------------------------------
void Simulation::execute()
{
  const unsigned int nGenerations = mParameters.nGenerations;
  mTimeSeries.timePoints.resize(nGenerations);
  for (unsigned int generation=0; generation<nGenerations; ++generation)
  {
    //log("Start of species selection. Nmales: " + IntToStr(mMales.size()) + ",Nfemales: " + IntToStr(mFemales.size()));
    mTimeSeries.timePoints[generation].speciesTallyOffspring.tallySpecies(mFemales,mMales);
    mTimeSeries.timePoints[generation].getSample(mFemales,mMales);

    speciesSelection();   //Selection on species, hybrids have disadvantage

    //log("Start of viabilitySelection. Nmales: " + IntToStr(mMales.size()) + ",Nfemales: " + IntToStr(mFemales.size()));
    mTimeSeries.timePoints[generation].speciesTallyAfterSpeciesSelection.tallySpecies(mFemales,mMales);

    viabilitySelection(); //Selection on trait, individuals with high trait have disadvantage

    //log("Selection before mating. Nmales: " + IntToStr(mMales.size()) + ",Nfemales: " + IntToStr(mFemales.size()));
    mTimeSeries.timePoints[generation].speciesTallyAfterTraitSelection.tallySpecies(mFemales,mMales);

    bool canDoSelection = mDensityDependentSelection->canDoSelection(mFemales,mMales,mParameters,mError);
    if (canDoSelection==false) break;
    mDensityDependentSelection->doSelection(mFemales,mMales,mParameters);

    mTimeSeries.timePoints[generation].speciesTallyAfterDensityDependentSelection.tallySpecies(mFemales,mMales);

    Offspring offspring = mMatingSystem->mate(mMales,mFemales,mParameters,mFemaleSampling,mTimeSeries.timePoints[generation].mateTally);

    std::swap(mFemales, offspring.females);
    std::swap(mMales, offspring.males);

    if (mFemales.size()==0) break;
    if (mMales.size()==0) break;

  }
}
//---------------------------------------------------------------------------
//Selection on species, hybrids have disadvantage
void Simulation::speciesSelection()
{
  //Easy implementation: hybrids just die
  for (unsigned int male = 0; male < mMales.size(); ++male)
  {
    while(male < mMales.size() && willDieSpecies(mMales[male])==true)
    {
      #ifdef NDEBUG
        mMales.erase(&mMales[male]);
      #else
        mMales.erase(mMales.begin() + male);
      #endif
    }
  }

  for (unsigned int female = 0; female < mFemales.size(); ++female)
  {
    while(female < mFemales.size() && willDieSpecies(mFemales[female])==true)
    {
      #ifdef NDEBUG
        mFemales.erase(&mFemales[female]);
      #else
        mFemales.erase(mFemales.begin() + female);
      #endif
    }
  }
}
//---------------------------------------------------------------------------
//Selection on trait, individuals with high trait have disadvantage
void Simulation::viabilitySelection()
{
  //Easy implementation: hybrids just die
  for (unsigned int male = 0; male < mMales.size(); ++male)
  {
    while(male < mMales.size() && willDieTrait(mMales[male])==true)
    {
      #ifdef NDEBUG
        mMales.erase(&mMales[male]);
      #else
        mMales.erase(mMales.begin() + male);
      #endif
    }
  }

  for (unsigned int female = 0; female < mFemales.size(); ++female)
  {
    while(female < mFemales.size() && willDiePreference(mFemales[female])==true)
    {
      #ifdef NDEBUG
        mFemales.erase(&mFemales[female]);
      #else
        mFemales.erase(mFemales.begin() + female);
      #endif
    }
  }
}
//---------------------------------------------------------------------------
//Kills the hybrids
bool Simulation::willDieSpecies(const Bird& bird)
{
  assert(bird.descent>=-1.0 && bird.descent<=1.0);
  const double chanceToSurvive
    = chanceToSurviveSpecies(bird.descent, mParameters.surviveSpeciesAlpha, mParameters.surviveSpeciesBeta);
  //const double chanceToSurvive = 1.0 - (mParameters.surviveSpeciesAlpha *
  //  std::exp(-mParameters.surviveSpeciesBeta * bird.descent * bird.descent));
  assert(chanceToSurvive>=0.0 && chanceToSurvive<=1.0);
  //Dot("bird.descent: " + FloatToStr(bird.descent)
  //  + ", chance to survive: " + FloatToStr(chanceToSurvive));
  if (rnd::uniform() > chanceToSurvive) return true;
  else return false;
}
//---------------------------------------------------------------------------
bool Simulation::willDieTrait(const Male& male)
{
  const double chanceToSurvive = chanceToSurviveTrait(male.trait, mParameters.costTrait);
  //const double chanceToSurvive
  //  = std::exp(-mParameters.costTrait * male.trait * male.trait);
  assert(chanceToSurvive>=0.0 && chanceToSurvive<=1.0);
  if (rnd::uniform() > chanceToSurvive) return true;
  else return false;
}
//---------------------------------------------------------------------------
bool Simulation::willDiePreference(const Female& female)
{
  const double chanceToSurvive
  //  = std::exp(-mParameters.costPreference * female.preference * female.preference);
    = chanceToSurvivePreference(female.preference, mParameters.costPreference);
  assert(chanceToSurvive>=0.0 && chanceToSurvive<=1.0);
  if (rnd::uniform() > chanceToSurvive) return true;
  else return false;
}
//---------------------------------------------------------------------------
void Simulation::setMatingSystem(
  const enumMatingSystem& matingSystem,
  const enumDensityDependentSelection& selection)
{
  setDensityDependentSelection(selection);

  switch(selection)
  {
    //Selection after mating
    case afterMating:
        switch(matingSystem)
        {
          case monogamy:
            mMatingSystem.reset(new StateMatingSystemMonogamyFixedNumberOffspring);
            break;
          case polygyny:
            mMatingSystem.reset(new StateMatingSystemPolygynyFixedNumberOffspring);
            break;
          default: assert(!"Unknown mating system"); std::exit(1);
        }
        break;
    //Selection before mating
    case beforeMating:
        switch(matingSystem)
        {
          case monogamy:
            mMatingSystem.reset(new StateMatingSystemMonogamyFreeNumberOffspring);
            break;
          case polygyny:
            mMatingSystem.reset(new StateMatingSystemPolygynyFreeNumberOffspring);
            break;
          default: assert(!"Unknown mating system"); std::exit(1);
        }
        break;
    default:
      assert(!"Unknown enumDensityDependentSelection"); std::exit(1);
  }
  Dot("Mating system set to: " + mMatingSystem->getMatingSystem());
}
//---------------------------------------------------------------------------
void Simulation::setFemaleSampling(const enumFemaleSampling& femaleSampling)
{
  switch(femaleSampling)
  {
    case bestOfNconspicific:
      mFemaleSampling.reset(new StateFemaleSamplingBestOfNconspicific);
      break;
    case bestOfNextremeTrait:
      mFemaleSampling.reset(new StateFemaleSamplingBestOfNextremeTrait);
      break;
    case bestOfNclosestTrait:
      mFemaleSampling.reset(new StateFemaleSamplingBestOfNclosestTrait);
      break;
    case fixedThresholdConspicific:
      mFemaleSampling.reset(new StateFemaleSamplingFixedThresholdConspicific);
      break;
    case fixedThresholdTraitSign:
      mFemaleSampling.reset(new StateFemaleSamplingFixedThresholdTraitSign);
      break;
    case fixedThresholdProbabilistic:
      mFemaleSampling.reset(new StateFemaleSamplingFixedThresholdProbabilistic);
      break;
    default:
      assert(!"Unknown female sampling"); std::exit(1);
  }
  Dot("Female sampling set to: " + mFemaleSampling->getFemaleSampling());
}
//---------------------------------------------------------------------------
void Simulation::setDensityDependentSelection(const enumDensityDependentSelection& selection)
{
  switch(selection)
  {
    case beforeMating:
      mDensityDependentSelection.reset(new StateDensityDependentSelectionBeforeMating);
      break;
    //case afterMating:
    //  mDensityDependentSelection.reset(new StateDensityDependentSelectionAfterMating);
    //  break;
    default:
      assert(!"Unknown enumDensityDependentSelection"); std::exit(1);
  }
  Dot("Density dependent selection set to: " + mDensityDependentSelection->getString());
}
//---------------------------------------------------------------------------
inline double Simulation::chanceToSurviveSpecies(const double& descent, const double& alpha, const double& beta)
{
  return 1.0 - (alpha * std::exp(-beta * descent * descent));
}
//---------------------------------------------------------------------------
inline double Simulation::chanceToSurviveTrait(const double& trait, const double& costTrait)
{
  return  std::exp(-costTrait * trait * trait);
}
//---------------------------------------------------------------------------
inline double Simulation::chanceToSurvivePreference(const double& preference, const double& costPreference)
{
  return  std::exp(-costPreference * preference * preference);
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

speciestally.h

 

#ifndef SPECIESTALLY_H
#define SPECIESTALLY_H
//---------------------------------------------------------------------------
#include <string>
#include <vector>
#include "bird.h"
//---------------------------------------------------------------------------
class SpeciesTally
{
  public:
  SpeciesTally() { reset(); }
  SpeciesTally(const std::vector<Female>& females, const std::vector<Male>& males)
  {
    reset();
    tallySpecies(females,males);
  }
  void reset();
  bool isNull() const;
  void tallySpecies(const std::vector<Female>& females, const std::vector<Male>& males);
  void operator+=(const SpeciesTally& speciesTally);
  void operator/=(const unsigned int& intValue);

  //Plain retrieval
  unsigned int getNmalesA()         const { return nMalesA;         }
  unsigned int getNmalesB()         const { return nMalesB;         }
  unsigned int getNfemalesA()       const { return nFemalesA;       }
  unsigned int getNfemalesB()       const { return nFemalesB;       }
  //Some group retrieval
  unsigned int getNallMalesA()   const { return nMalesA; }
  unsigned int getNallMalesB()   const { return nMalesB; }
  unsigned int getNallMales()    const { return nMalesA + nMalesB; }
  unsigned int getNallFemalesA() const { return nFemalesA; }
  unsigned int getNallFemalesB() const { return nFemalesB; }
  unsigned int getNallFemales()  const { return nFemalesA + nFemalesB; }
  unsigned int getNall()         const { return getNallFemales() + getNallMales(); }

  private:
  //The private variables it is all about
  unsigned int nMalesA;         //SpeciesValue < -0.5
  unsigned int nMalesB;         //SpeciesValue > 0.5
  unsigned int nFemalesA;       //SpeciesValue < -0.5
  unsigned int nFemalesB;       //SpeciesValue > 0.5
};
//---------------------------------------------------------------------------


#endif // SPECIESTALLY_H

 

 

 

 

 

speciestally.cpp

 

//---------------------------------------------------------------------------
#include "speciestally.h"
//---------------------------------------------------------------------------
void SpeciesTally::tallySpecies(const std::vector<Female>& females, const std::vector<Male>& males)
{
  //Females first
  const unsigned int nFemales = females.size();
  for (unsigned i=0; i<nFemales; ++i)
  {
    const enumSpecies species = females[i].species;
    if (species==piedFlycatcher) ++(nFemalesA);
    else                         ++(nFemalesB);
  }
  //Tally the men
  const unsigned int nMales = males.size();
  for (unsigned i=0; i<nMales; ++i)
  {
    const enumSpecies species = males[i].species;
    if (species==piedFlycatcher) ++(nMalesA);
    else                         ++(nMalesB);
  }

  /*
  //Females first
  const unsigned int nFemales = females.size();
  for (unsigned i=0; i<nFemales; ++i)
  {
    const double descent = females[i].descent;
    if (descent < -0.5)     ++(nFemalesA);
    else if (descent < 0.0) ++(nFemaleHybridsA);
    else if (descent < 0.5) ++(nFemaleHybridsB);
    else                         ++(nFemalesB);
  }
  //Tally the men
  const unsigned int nMales = males.size();
  for (unsigned i=0; i<nMales; ++i)
  {
    const double descent = males[i].descent;
    if (descent < -0.5)     ++(nMalesA);
    else if (descent < 0.0) ++(nMaleHybridsA);
    else if (descent < 0.5) ++(nMaleHybridsB);
    else                         ++(nMalesB);
  }
  */
}
//---------------------------------------------------------------------------
void SpeciesTally::reset()
{
  nMalesA = 0;
  nMalesB = 0;
  nFemalesA = 0;
  nFemalesB = 0;
}
//---------------------------------------------------------------------------
bool SpeciesTally::isNull() const
{
  if (nMalesA>0) return false;
  if (nMalesB>0) return false;
  if (nFemalesA>0) return false;
  if (nFemalesB>0) return false;
  return true;
}
//---------------------------------------------------------------------------
void SpeciesTally::operator+=(const SpeciesTally& speciesTally)
{
  nMalesA+=speciesTally.nMalesA;
  nMalesB+=speciesTally.nMalesB;
  nFemalesA+=speciesTally.nFemalesA;
  nFemalesB+=speciesTally.nFemalesB;
}
//---------------------------------------------------------------------------
void SpeciesTally::operator/=(const unsigned int& intValue)
{
  nMalesA/=intValue;
  nMalesB/=intValue;
  nFemalesA/=intValue;
  nFemalesB/=intValue;
}
//---------------------------------------------------------------------------

 

 

 

 

 

timepoint.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef TIMEPOINT_H
#define TIMEPOINT_H
//---------------------------------------------------------------------------
#include <vector>
#include <string>
#include <algorithm>
#include "matetally.h"
#include "speciestally.h"
//---------------------------------------------------------------------------
class TimePoint
{
  public:
  TimePoint() { reset(); }
  std::vector<double> traits;
  std::vector<double> preferences;
  std::vector<double> descents;
  MateTally mateTally;
  SpeciesTally speciesTallyOffspring;
  SpeciesTally speciesTallyAfterSpeciesSelection;
  SpeciesTally speciesTallyAfterTraitSelection;
  SpeciesTally speciesTallyAfterDensityDependentSelection;

  void reset();
  bool isNull() const;
  void getSample(const std::vector<Female>& females, const std::vector<Male>& males);

  void operator/=(const unsigned int& valueInt) //For taking the average of multiple timepoints
  {
    //const double value = static_cast<double>(valueInt);
    speciesTallyOffspring/=valueInt;
    speciesTallyAfterSpeciesSelection/=valueInt;
    speciesTallyAfterTraitSelection/=valueInt;
    speciesTallyAfterDensityDependentSelection/=valueInt;
    std::random_shuffle(traits.begin(),traits.end());
    std::random_shuffle(preferences.begin(),preferences.end());
    std::random_shuffle(descents.begin(),descents.end());
    traits.resize(traits.size()/valueInt);
    preferences.resize(preferences.size()/valueInt);
    descents.resize(descents.size()/valueInt);

  }
  void operator+=(const TimePoint& timePoint) //For taking the average of multiple timepoints
  {
    mateTally+=timePoint.mateTally;
    speciesTallyOffspring+=timePoint.speciesTallyOffspring;
    speciesTallyAfterSpeciesSelection+=timePoint.speciesTallyAfterSpeciesSelection;
    speciesTallyAfterTraitSelection+=timePoint.speciesTallyAfterTraitSelection;
    speciesTallyAfterDensityDependentSelection+=timePoint.speciesTallyAfterDensityDependentSelection;

    const unsigned int traitSize = timePoint.traits.size();
    for (unsigned int i=0; i<traitSize; ++i) traits.push_back(timePoint.traits[i]);
    const unsigned int preferenceSize = timePoint.preferences.size();
    for (unsigned int i=0; i<preferenceSize; ++i) preferences.push_back(timePoint.preferences[i]);
    const unsigned int descentSize = timePoint.descents.size();
    for (unsigned int i=0; i<descentSize; ++i) descents.push_back(timePoint.descents[i]);
  }
};
//---------------------------------------------------------------------------
#endif // TIMEPOINT_H

 

 

 

 

 

timepoint.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "timepoint.h"
//---------------------------------------------------------------------------
void TimePoint::reset()
{
  traits.resize(0);
  preferences.resize(0);
  descents.resize(0);
  mateTally.reset();
  speciesTallyOffspring.reset();
  speciesTallyAfterSpeciesSelection.reset();
  speciesTallyAfterTraitSelection.reset();
  speciesTallyAfterDensityDependentSelection.reset();
}
//---------------------------------------------------------------------------
bool TimePoint::isNull() const
{
  if (traits.size()>0) return false;
  if (preferences.size()>0) return false;
  if (descents.size()>0) return false;
  if (mateTally.isNull()==false) return false;
  if (speciesTallyOffspring.isNull()==false) return false;
  if (speciesTallyAfterSpeciesSelection.isNull()==false) return false;
  if (speciesTallyAfterTraitSelection.isNull()==false) return false;
  if (speciesTallyAfterDensityDependentSelection.isNull()==false) return false;
  return true;
}
//---------------------------------------------------------------------------
void TimePoint::getSample(const std::vector<Female>& females, const std::vector<Male>& males)
{
  const unsigned int sampleSize = 2;
  //Get sample of female preferences
  unsigned int nFemalesA = 0, nFemalesB = 0;
  const unsigned int nFemales = females.size();
  for (unsigned int i=0; i<nFemales ; ++i)
  {
    if (females[i].species==piedFlycatcher)
    {
      if (nFemalesA<sampleSize)
      {
        preferences.push_back(females[i].preference);
        descents.push_back(females[i].descent);
        ++nFemalesA;
      }
    }
    else
    {
      if (nFemalesB<sampleSize)
      {
        preferences.push_back(females[i].preference);
        descents.push_back(females[i].descent);
        ++nFemalesB;
      }
    }
    if (nFemalesA==sampleSize && nFemalesB==sampleSize) break;
    assert(nFemalesA<=sampleSize);
    assert(nFemalesB<=sampleSize);
  }

  //Get sample of male traits
  unsigned int nMalesA = 0, nMalesB = 0;
  const unsigned int nMales = males.size();
  for (unsigned int i=0; i<nMales ; ++i)
  {
    if (males[i].species==piedFlycatcher)
    {
      if (nMalesA<sampleSize)
      {
        traits.push_back(males[i].trait);
        descents.push_back(males[i].descent);
        ++nMalesA;
      }
    }
    else
    {
      if (nMalesB<sampleSize)
      {
        traits.push_back(males[i].trait);
        descents.push_back(males[i].descent);
        ++nMalesB;
      }
    }
    if (nMalesA==sampleSize && nMalesB==sampleSize) break;
    assert(nMalesA<=sampleSize);
    assert(nMalesB<=sampleSize);
  }
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

 

 

 

 

 

timeseries.h

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#ifndef TIMESERIES_H
#define TIMESERIES_H
//---------------------------------------------------------------------------
#include <vector>
#include <assert.h>
#include "timepoint.h"
//---------------------------------------------------------------------------
class TimeSeries
{
  public:
  std::vector<TimePoint> timePoints;
  bool isNull() const
  {
    if (timePoints.size()>0) return false;
    return true;
  }

  void operator+=(const TimeSeries& timeSeries) //For taking the average of multiple timepoints
  {
    const unsigned int nTimePoints = timePoints.size();
    assert(timeSeries.timePoints.size()==nTimePoints);
    for (unsigned i=0; i<nTimePoints; ++i)
    {
      timePoints[i]+=timeSeries.timePoints[i];
    }
  }
  void operator/=(const unsigned int& valueInt) //For taking the average of multiple timepoints
  {
    const unsigned int nTimePoints = timePoints.size();
    for (unsigned i=0; i<nTimePoints; ++i)
    {
      timePoints[i]/=valueInt;
    }
  }
};
//---------------------------------------------------------------------------
void getMeanAndStdErrorEndPoint(const std::vector<TimeSeries>& timeSeries, TimePoint& mean, TimePoint& stdError);
void getMeanAndStdErrorTimeSeries(const std::vector<TimeSeries>& timeSeries, TimeSeries& mean, TimeSeries& stdError);
//TimeSeries getAverageTimeSeries(const std::vector<TimeSeries>& timeSeries);
//TimePoint getAverageEndTimePoint(const std::vector<TimeSeries>& timeSeries);
//---------------------------------------------------------------------------
#endif // TIMESERIES_H

 

 

 

 

 

timeseries.cpp

 

//---------------------------------------------------------------------------
/*
  The Chris Wiley Project, simulation on mixed-pair mating
  Copyright (C) 2009  Richl Bilderbeek

  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/>.
*/
//---------------------------------------------------------------------------
// From http://www.richelbilderbeek.nl
//---------------------------------------------------------------------------
#include "timeseries.h"
//---------------------------------------------------------------------------
void getMeanAndStdErrorEndPoint(const std::vector<TimeSeries>& timeSeries, TimePoint& mean, TimePoint& stdError)
{
  assert(mean.isNull()==true);
  assert(stdError.isNull()==true);

  //Only set: biasA, biasB, fractionMixedPairs

  const unsigned int nTimeSeries = timeSeries.size();
  for (unsigned i=0; i<nTimeSeries; ++i)
  {
    //Get the index of the final index
    const unsigned int time = timeSeries[i].timePoints.size() - 1;
    mean += timeSeries[i].timePoints[time];
  }
  mean/=nTimeSeries;
}
//---------------------------------------------------------------------------
void getMeanAndStdErrorTimeSeries(
  const std::vector<TimeSeries>& timeSeries,
  TimeSeries& mean,
  TimeSeries& /*stdError*/)
{
  assert(mean.isNull()==true);
  mean = timeSeries[0];
  const unsigned int nTimeSeries = timeSeries.size();
  for (unsigned i=1; i<nTimeSeries; ++i)
  {
    //Get the index of the final index
    mean += timeSeries[i];
  }
  mean/=nTimeSeries;
}
//---------------------------------------------------------------------------

 

 

 

 

 

Go back to Richel Bilderbeek's C++ page.

Go back to Richel Bilderbeek's homepage.

 

Valid XHTML 1.0 Strict

This page has been created by the tool CodeToHtml