README for SimpleCrowdingModelFlatTextArchive
This is the source code and supporting data files for the "Simple Crowding Model" used to generate the fits 
in Kieras(2020) technical report:

Kieras, D. (2020). Visual Search Explained with a Computational Cognitive Architecture: Early
Visual Processes, Eye Movements, and Task Strategies. 
Technical Report TR-20/ONR-EPIC-21, 
University of Michigan Electrical Engineering and Computer Science Department.
August 20, 2020. DTIC AD1114079 DOI: 10.13140/RG.2.2.23716.91520.

The code was written in Standard C++ 17. It should need at most only minor updating to run in newer
versions of C++. Separate out the individual files and build the executable using an IDE, a simple makefile or even simply compile *.cpp to produce the executable.

The top level main() module is CrowdingSim.cpp. 

The major submodules are Program_constants, Display_creator, Search_object, Visual_model, Strategy, 
PrdObsStatistics, and Trial_statistics. 

The model is configured for a particular condition, strategy, and visual parameters by editing the contents
of Program_constants.cpp.

The program reads the observed data values from one of the files in RunTimeInputFiles, and uses these
to compute goodness of fit statistics and output a table that can be copy-pasted into Excel to produce
graphs of the predicted and observed results.

The model code here has copies of components from the EPIC Architecture Simulation System.
This was used in order to maximize code sharing between this specialized 
but efficient model-fitting software and the general production-system-based EPIC architecture. 
This is the source of the many .h,.cpp file pairs such as Geometry, Random_utilities, etc.

The observed input data files are those with names starting with "ObsData". These data were 
derived from the data files made available by Wolfe at 
Visual Attention Lab (n.d.). Data Sets and look up Stimulus Sets: What can 8000 trials tell you
about visual search? Retrieved April 23, 2018 from http://search.bwh.harvard.edu/new/data_set_files.html.

------------------------------------------------------------
// Assert_throw.h
#ifndef ASSERT_THROW_H
#define ASSERT_THROW_H

#include "Exception.h"
#include <string>

// this class always declared
class Assertion_exception : public Exception {
public:
	Assertion_exception (const char * condition_, const char * filename_, int lineno_);
	~Assertion_exception() throw()
		{}
	const char * what() const throw() override
		{return msg.c_str();}
private:
	std::string msg;
};

// turn macro on or off depending on NDEBUG

#ifdef NDEBUG
#define Assert(condition) ((void) 0)
#else
#define Assert(condition) ((condition) ? ((void) 0) : throw Assertion_exception(#condition,__FILE__, __LINE__))
#endif

#endif

------------------------------------------------------------
// Assert_throw.cpp
#include "Assert_throw.h"

#include <string>
#include <sstream>

using std::string;
using std::ostringstream;

Assertion_exception::Assertion_exception (const char * condition_, const char * filename_, int lineno_)
{
	ostringstream ss;
	ss << "Assertion exception: " << condition_ << " failed in " << filename_ << " line " << lineno_;
	msg = ss.str();
}

------------------------------------------------------------
// CrowdingSimMain.cpp
//  
// Simple Crowding Model 
//  Created by David Kieras on 3/7/18.
//

// as of 3/9/18, reproduced number of fixations produced by EPIC model of same date (used in ICCM2018Sub)

#include "Program_constants.h"
#include "Search_object.h"
//#include "Visual_modelV2b.h"
#include "Visual_modelV2e.h"
#include "Display_creator.h"
#include "Strategy.h"
#include "Trial_statistics.h"
#include "PrdObsStatistics.h"

/*#include "EPICLib/Symbol.h"
#include "EPICLib/Geometry.h"
#include "EPICLib/Statistics.h"
#include "EPICLib/Ostream_format_saver.h"
#include "EPICLib/Random_utilities.h"
#include "EPICLib/Symbol.h" */
#include "Geometry.h"
#include "Statistics.h"
#include "Ostream_format_saver.h"
#include "Random_utilities.h"

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <array>
#include <vector>
#include <algorithm>
#include <cassert>
#include <cfloat>
#include <limits>
#include <exception>


namespace GU = Geometry_Utilities;

using namespace std;

void do_grid_search();
void run_experiment();
void do_condition(Display_creator& display_creator, Condition_e cond, Polarity_e polarity, int set_size, int n_trials, PrdObsStatistics& prdobsstatistics);
//void do_trial(Display_creator& display_creator, Condition_e cond, int set_size, Polarity_e polarity, Statistics_accumulator& statistics);
//GU::Point move_eye(GU::Point current_eye_location, GU::Point obj_location);
//bool make_positive_response(Polarity_e polarity, Statistics_accumulator& statistics);
//bool make_negative_response(Polarity_e polarity, Statistics_accumulator& statistics);

const array<const char* const, 3> cond_labels{"CSF", "COC", "SHP"};
const array<const char* const, 2> pol_labels{"Neg", "Pos"};

// experiment mode runs conditions and writes predicted values and parameter+gof values to cout.
// grid-search mode runs conditions and writes experiment-mode output to a file and parameter+gof values to cout.

// try as global
//PrdObsStatistics prdobsstatistics;
ofstream output_file;

// following declared in Program_constants.h


//bool cout_experiment_results = true; // set to false for compact grid search results

int main()
{
    Ostream_format_saver format_save(cout);
    // Set output to show two decimal places for floating point numbers
//    cout << fixed << setprecision(2) << endl;
    cout.setf(ios::fixed, ios::floatfield);
    cout.precision(3);
    
/*    cout << "Normal cdf test mean = 0, sd = 1" << endl;
    for(double x = -6; x < 6.; x = x+1.0)
        cout << x << ' ' << normal_cdf(x, 0.0, 1.0) << endl;
 
    cout << "Normal cdf test for mean = 10, sd = 5" << endl;
    for(double x = -20; x < 20.; x = x+1.0)
        cout << x << ' ' << normal_cdf(x, 10., 5.0) << endl;
 */
/*    cout << "Gaussian detection probability test with a = 0, b = .5, sd = 0.5, size = GU::Size(3.5, 1.)" << endl;
    for (double ecc = 0.; ecc < 20.; ecc = ecc + 1)
        cout << ecc << ' ' << LinearConstSD_detection_probability(ecc, GU::Size(3.5, 1.), 0.0, 0.5, 0.5) << endl;
*/

//    cout << "_DECIMAL DIG: double precision: " << DBL_DECIMAL_DIG << " long double precision: " << LDBL_DECIMAL_DIG << endl;
//    cout << "sizeof double: " << sizeof(double) << " sizeof long double: " << sizeof(long double) << endl;
 /*   // load the desired data file
    ifstream infile(infile_name_c);
    if(!infile.good()) {
        cout << "Could not open observed data file \"" << infile_name_c << "\"!" << endl;
        return 1;
        }
    cout << "Observed data from " << infile_name_c << "\n";
    prdobsstatistics.load_obs_data(infile, false); // second parameter true means echo_print
*/
    cout << "Observed data from " << infile_name_c << "\n";
    cout << n_trials_c << " trials/condition" << endl;
    cout << Strategy::get_strategy_id_string();
    cout << Visual_model::get_id_string() << endl;
//    cout << ((apply_crowding_effect) ? "Crowding applied" : "No Crowding effect") << endl;

    try {
    switch(run_mode_c) {
        case Run_mode_e::EXPERIMENT:
            run_experiment();
            break;
        case Run_mode_e::GRID_SEARCH: {
            output_file.open("Outputfile.txt");
            if(!output_file.is_open())
                throw runtime_error("Could not open \"Outputfile.txt\"");
            do_grid_search();
            break;
            }
        default:
            assert(!"Invalid run_mode");
            break;
        }
/*    if(run_experiment_c)
        run_experiment();
    else
        do_grid_search();
*/
    }
    catch(exception& x) {
        cout << "Exception: " << x.what() << endl;
        }
}

const char* const model_output_headings_c = "\nCnd\t\tPol\tSS\tCN\tCRT\tEN\tERT\tER\tCMFix\tEMFix\tMITF\tMRptF\tMTgtF\tMNTID\tMNTIB\tMNIT";
const char* const parameter_metrics_heading_c =
    "C_av\tC_cp\tO_av\tO_cp\tS_av\tS_cs\tS_cp\tVDly\tPbER\tNbER\tMaxNfix\tCPskip\tCNskip\tMaxITF"
    "\tRT: Rsq\taae\taare\tFoM\tER: Rsq\taae\taare\tFoMa\tFoMr\tWAFoMs\tN";


void do_grid_search()
{
    cout << parameter_metrics_heading_c << endl;
    if(output_file.is_open())
        output_file << parameter_metrics_heading_c << endl;

//    cout << "C_av\tC_cp\tO_av\tO_cp\tS_av\tS_cs\tS_cp\tVDly\tPbER\tNbER\tMaxNfix\t";
//    cout << "RT: Rsq\taae\taare\tFoM\tER: Rsq\taae\taare\tFoMa\tFoMr\tWAFoMs\tN\n";
//    for(max_n_fixations_g = 2; max_n_fixations_g <= 3; max_n_fixations_g += 1) {
        for(color_slope_g = 0.0; color_slope_g <= 0.4; color_slope_g += 0.05) {
//            for(color_crowding_probability_g = .1; color_crowding_probability_g <= 0.8; color_crowding_probability_g += 0.1) {
//            for(orientation_slope_g = 0.05; orientation_slope_g <= 0.20; orientation_slope_g += 0.05) {
//            for(orientation_crowding_probability_g = .0; orientation_crowding_probability_g <= .1; orientation_crowding_probability_g += 0.05) {
 //                   for(shape_slope_g = .3; shape_slope_g <= .4; shape_slope_g += 0.025) {
 //                   for(shape_crowding_probability_g = 0.00; shape_crowding_probability_g <= .05; shape_crowding_probability_g += 0.005) {
                            run_experiment();
//                        } // shape_crowding_probability_g
//                        } // shape_slope_g
//                } //orientation_crowding_probability_g
//                } // orientation_slope_g
//            } // color_crowding_probability_g
            } // color_slope_g
//        } // max_n_fixations_g
}


void run_experiment()
{
    PrdObsStatistics prdobsstatistics;
    ifstream infile(infile_name_c);
    if(!infile.good()) {
        throw runtime_error("Could not open observed data file!");
//        cout << "Could not open observed data file \"" << infile_name_c << "\"!" << endl;
//        return 1;
        }
    prdobsstatistics.load_obs_data(infile, false); // second parameter true means echo_print
//    prdobsstatistics.reset();
    reset_random_number_generator_seed();
    Display_creator display_creator;

    if (run_mode_c == Run_mode_e::EXPERIMENT) {
        cout << model_output_headings_c << endl;
        }
    if(output_file.is_open()) {
        output_file << model_output_headings_c << endl;
        }


//    const auto conditions = {Condition_e::CSF, Condition_e::COC, Condition_e::SHP};
//    const auto conditions = {Condition_e::CSF, Condition_e::COC};
//    const auto conditions = {Condition_e::CSF};
//    const auto conditions = {Condition_e::COC};
//    const auto conditions = {Condition_e::SHP};

//    const auto polarities = {Polarity_e::NEGATIVE, Polarity_e::POSITIVE};
//    const auto polarities = {Polarity_e::POSITIVE};

//    const auto set_sizes = {3,6,12,18};
//    const auto set_sizes = {12};

    for(auto cond : conditions) {
        for(auto polarity : polarities) {
          for(int set_size : set_sizes) {
                do_condition(display_creator, cond, polarity, set_size, n_trials_c, prdobsstatistics);
                }
            }
        }

    ostringstream oss;
    oss << fixed << setprecision(3);
//    oss << parameter_metrics_heading_c << endl;
    oss << color_slope_g << "\t" << color_crowding_probability_g << "\t";
    oss << orientation_slope_g << "\t" << orientation_crowding_probability_g << "\t";
//    cout << shape_slope_g << "\t" << shape_crowding_probability_g << "\t";
    oss << shape_slope_g << "\t" << shape_crowding_slope_g << "\t" << shape_crowding_probability_g << "\t";
    oss << visual_delay_time_g << "\t";
    oss << setprecision(4) << positive_baseline_error_rate_c << "\t" << negative_baseline_error_rate_c << "\t";
    oss << setprecision(3);
    // if mixture probability > 0.0 AND max_n_fixations < 99, then use _g value,
    // otherwise calculate average value using mixture probability
    //
    if(max_n_fixations_mixture_probability_c > 0.0 && max_n_fixations_g < 99)
        oss << calculate_average_max_n_fixations() << "\t";
    else
        oss << max_n_fixations_g << "\t";;
    oss << close_enough_to_bypass_confirm_positive_c << "\t";
    oss << enough_fixations_to_bypass_confirm_negative_c << "\t";
    oss << too_many_illusory_targets_c << "\t";
    prdobsstatistics.output(oss);
    
     if (run_mode_c == Run_mode_e::EXPERIMENT) {
        cout << parameter_metrics_heading_c << endl;
        }
   // always output parameters and metrics
    cout << oss.str() << endl;
    if(output_file.is_open())
        output_file << oss.str() << endl;
}


void do_condition(Display_creator& display_creator, Condition_e cond, Polarity_e polarity, int set_size, int n_trials, PrdObsStatistics& prdobsstatistics)
{
    Trial_statistics statistics;
    Strategy strategy(display_creator, cond, set_size, polarity, statistics);
	
    display_creator.reset_display_statistics();
    
    // trial needs to be a global variable so that debugging output can include it
 
    for(trial_g = 0; trial_g < n_trials; trial_g++) {
        if(strategy_trace_output) cout << "\nTrial " << trial_g << endl;
        strategy.do_trial();
        }
    ostringstream oss;
    oss << fixed << setprecision(3);
    oss  << cond_labels[static_cast<int>(cond)] << "\t\t" << pol_labels[static_cast<int>(polarity)] << '\t' << set_size << '\t'
            << statistics.mean_RT.get_n() << '\t'
            << statistics.mean_RT.get_mean() << '\t'
            << statistics.mean_error_RT.get_n() << '\t'
            << statistics.mean_error_RT.get_mean() << '\t'
            << statistics.prop_errors.get_proportion() << '\t'
            << statistics.mean_n_fixations.get_mean() << '\t'
            << statistics.mean_error_n_fixations.get_mean() << '\t'
            << statistics.mean_illusory_target_fixations.get_mean() << '\t'
            << statistics.mean_repeat_fixations.get_mean() << '\t'
            << statistics.mean_target_fixations.get_mean() << '\t'
            << statistics.mean_n_target_illusory_distractor.get_mean() << '\t'
            << statistics.mean_n_target_illusory_blank.get_mean() << '\t'
            << statistics.mean_n_illusory_targets.get_mean()
//            << statistics.mean_n_target_sens_property_present.get_mean() << '\t'
//            << statistics.mean_n_target_perc_property_present.get_mean() << '\t'
//            << statistics.mean_n_target_property_overwritten.get_mean()
            ;

    if (run_mode_c == Run_mode_e::EXPERIMENT) {
        cout << oss.str() << endl;
/*    	cout << cond_labels[static_cast<int>(cond)] << "\t\t" << pol_labels[static_cast<int>(polarity)] << '\t' << set_size << '\t'
        << statistics.mean_RT.get_n() << '\t'
        << statistics.mean_RT.get_mean() << '\t'
        << statistics.prop_errors.get_proportion() << '\t'
        << statistics.mean_n_fixations.get_mean() << '\t'
        << statistics.mean_illusory_target_fixations.get_mean() << endl;
*/
        }
    if(output_file.is_open()) {
        output_file << oss.str() << endl;
        }
    if(compute_and_output_display_statistics)
        display_creator.output_display_statistics();

    prdobsstatistics.update(cond, polarity, set_size, statistics.mean_RT.get_mean(),statistics.prop_errors.get_proportion());
}



------------------------------------------------------------
// Display_creator.h
//  
//
//  Created by David Kieras on 3/7/18.
//

#ifndef DISPLAY_CREATOR_H
#define DISPLAY_CREATOR_H

#include "Search_object.h"
#include "Statistics.h"
//#include "EPICLib/Statistics.h"

#include <vector>

// a class for accumulating display characteristics relevant to crowding
struct Display_properties_accumulator {
    int n_trials;
    Mean_accumulator number_of_objects;	// for validity checking
    Mean_accumulator eccentricities;
    Mean_accumulator spacings;
    Mean_accumulator numbers_of_crowding_flankers;
    Mean_accumulator spacings_of_crowding_flankers;
    static constexpr int n_crowding_flankers_dbn_bins_c = 10;
    static constexpr int eccentricity_n_bins_c = 32;
    static constexpr int eccentricity_bin_size_c = 1;
    Distribution_accumulator dbn_of_eccentricities{eccentricity_n_bins_c, eccentricity_bin_size_c};
    // 10 bins for zero through 9 flankers, more aren't counted in last bin but included in total number
    Discrete_distribution_accumulator dbn_of_number_of_crowding_flankers{n_crowding_flankers_dbn_bins_c};
    // 32 bins of 1 deg each for the average number of crowders at that eccentricity
    Binned_mean_accumulators mean_n_crowders_for_eccentricity{eccentricity_n_bins_c, eccentricity_bin_size_c};
    void reset()
        {
        n_trials = 0;
        number_of_objects.reset();
        eccentricities.reset();
        spacings.reset();
        numbers_of_crowding_flankers.reset();
        spacings_of_crowding_flankers.reset();
        dbn_of_eccentricities.reset();
        dbn_of_number_of_crowding_flankers.reset();
        mean_n_crowders_for_eccentricity.reset();
        }
    // update using means and counts from another accumulator of the same type
    void update(const Display_properties_accumulator& source)
    {
        eccentricities.update(source.eccentricities.get_mean());
        spacings.update(source.spacings.get_mean());
        numbers_of_crowding_flankers.update(source.numbers_of_crowding_flankers.get_mean());
        spacings_of_crowding_flankers.update(source.spacings_of_crowding_flankers.get_mean());
        dbn_of_eccentricities.update(source.dbn_of_eccentricities);
        dbn_of_number_of_crowding_flankers.update(source.dbn_of_number_of_crowding_flankers);
        mean_n_crowders_for_eccentricity.update_with_bin_means(source.mean_n_crowders_for_eccentricity);
    }
    
    void output() const;
};


class Display_creator {
public:

    Display_creator();
    void initialize();
     // return a vector of Search_object, each fully constructed, positioned, and ready to be used
    Search_objects_t create_display(Condition_e condition, int set_size, Polarity_e polarity);
	void reset_display_statistics();	// call to reset accumulators (e.g. before each condition)
    void output_display_statistics() const;


private:
    // stimulus generation and representation
    using Locations_t = std::vector<GU::Point> ;
    Locations_t all_locations;
    

    Search_objects_t search_objects;
    
    Condition_e condition;
    Polarity_e polarity;
    int n_objects;
    
    Search_object target;
    
    
    // stimulus characteristic data calculated for each display
    // following is for characteristics measured from the fixation point -
    // corresponding to information available at the trial start, before any eye movements
    
    
    Display_properties_accumulator initial_display_properties;
    Display_properties_accumulator whole_display_properties;

    
    // helpers
    void generate_stimulus();
    void generate_CSF_stimulus();
    void generate_COC_stimulus();
    void generate_SHP_stimulus();
    void present_stimulus();
    void present_search_objects();
    void calculate_display_metrics();
    void calculate_display_statistics();
    void calculate_possible_next_fixation_statistics(GU::Point current_assumed_fixation_point, int possible_next_fixation_index, Display_properties_accumulator& properties_accumulator);
    void output_display_properties_accumulator(const char* message,const Display_properties_accumulator& properties_accumulator) const;
};


#endif /* DISPLAY_CREATOR_H */

------------------------------------------------------------
// Display_creator.cpp
//  
//
//  Created by David Kieras on 3/7/18.
//

/*
Wolfe, Palmer, Horowitz 2010 model
Conditions:
Number of objects: 3, 6, 12, 18
Cue condition: Color Single Feature (SCF), Color Orientation Conjunction (COC), Shape (SHP)
Trial Polarity: Positive (Cued target is present), Negative (cued target is absent).
50% Positive, 50% Negative

Stimuli:
On each trial, positive/negative chosen at random, number of objects chosen at random

CSF: vertical bars 1 X 3.5; target: red vertical bar; distractors: green vertical bar

COC: vertical 1 X 3.5 or horizontal 3.5 X 1 bars; target: red vertical; distractors: red horizontal, green vertical.

SHP: "digital" 2s and 5s 1.5 X 2.7; target: 2, distractors 5

Visual field is about 22.5 x 22.5 degrees, divided up into 4.5 X 4.5 degree squares
in 5 X 5 matrix.
Apparently objects were put into a randomly selected square with no more than one object is placed in a square.
An object within a square was put at a random location within the square (not further described).
This is location 'jitter'.

Say we will maintain at least 0.25 dva between objects then: gives 4 dva range both x and y center location to vary in

CSF: center +/- 1.5 x, +/- .25 y U(-1.5, + 1.5) + x, U(-.25, +.25) + y or
                                    uniform_random_variable(x, 1.5), (y, .25)
COC vertical bar, +/- 1.5 x, +/- .25 y
    horizontal bar +/- .25 x, +/- 1.5 y
SHP center +/- 1.25 x, .65 y

Trial sequence:

1. Central fixation point .7 X .7 dva. Subjects were told to keep eyes on fixation point
but eye movements were not monitored. Fixation point apparently remains - not clear

2. Beep, then 500 ms delay, display appears.

3. S presses 'z' or '/' for whether target is absent or present.

4. Display disappears.

5. Correct/error feedback fro 500 ms.

6. Wait ITI=1000 ms for next trial start.

Subject could pause experiment by pressing space bar

Stimulus construction:
Randomize locations, then put N distractors in;
flip coin for trial polarity;
if positive, replace first distractor with target (chosen from cue condition);
if negative, leave as is.

Data:
RT in each condition, count of number of eye movements after space bar
Dwell is actually not meaningful for first fixation (eye was placed at fixation point long ago)
or last fixation (response is underway, eye not moved until trial is over).
So fixations should not be counted until trial starts, and dwell statistics on them not computed until trial ends.


*/

#include "Display_creator.h"
#include "Search_object.h"
#include "Program_constants.h"
/*#include "EPICLib/Geometry.h"
#include "EPICLib/Numeric_utilities.h"
#include "EPICLib/Random_utilities.h"
#include "EPICLib/Symbol_utilities.h"
#include "EPICLib/Ostream_format_saver.h" */
#include "Geometry.h"
#include "Numeric_utilities.h"
#include "Random_utilities.h"
#include "Symbol_utilities.h"
#include "Ostream_format_saver.h"
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <algorithm>
#include <random>
#include <cassert>
#include <cmath>

namespace GU = Geometry_Utilities;
//using GU::Point;
using namespace std;

// handy internal constants
const Symbol Fixation_point_c("Fixation_point");
const Symbol Display_c("Display");
const Symbol Dummy_c("Dummy");
const Symbol Go_response_c("SP");
const Symbol Response_absent_c("z");
const Symbol Response_present_c("/");


const char Red_c{'R'};
const char Green_c{'G'};
const char Black_c{'B'};
const char Vertical_c{'V'};
const char Horizontal_c{'H'};
const char Rectangle_c{'r'};
const char Digital_2_c{'2'};
const char Digital_5_c{'5'};

const double CSF_x_jitter = 1.5;
const double CSF_y_jitter = .25;
const double SHP_x_jitter = 1.25;
const double SHP_y_jitter = .65;
const double COCH_x_jitter = .25;
const double COCH_y_jitter = 1.5;
const double COCV_x_jitter = 1.5;
const double COCV_y_jitter = .25;

GU::Point jitter(GU::Point point, double x_deviation, double y_deviation);


Display_creator::Display_creator()
{
//Visual field is about 22.5 x 22.5 degrees, divided up into 4.5 X 4.5 degree squares
//in 5 X 5 matrix.
    // fill the location vector 5 X 5 grid, skip the center position
    // locations are centers of 4.5 x 4.5 boxes; skip the center position
    // maybe draw a box?
    for(int x = -2; x <= +2; x++) {
        for(int y = -2; y <= +2; y++) {
            if(x == 0 && y == 0)
                continue;
            all_locations.push_back(GU::Point(x * 4.5, y * 4.5));
            }
        }
    
    initialize();

}



// called externally to start up device
void Display_creator::initialize()
{
    // reset accumulators
    whole_display_properties.reset();
    initial_display_properties.reset();

}


std::vector<Search_object> Display_creator::create_display(Condition_e condition_, int set_size, Polarity_e polarity_)
{
    condition = condition_;
    n_objects = set_size; // to reuse code and improve modularity
    polarity = polarity_;
    
    //Was done>TO DO: replace with modern std::shuffle
    // get a randomization of the possible locations
//    random_shuffle(all_locations.begin(), all_locations.end(), random_int);
    shuffle(all_locations.begin(), all_locations.end(), get_Random_engine());

    // empty the container of search objects
    search_objects.clear();
    
    switch (condition) {
        case Condition_e::CSF: generate_CSF_stimulus(); break;
        case Condition_e::COC: generate_COC_stimulus(); break;
        case Condition_e::SHP: generate_SHP_stimulus(); break;
        };
    // add any essential metrics to the search objects
    calculate_display_metrics();
    // calculate and average the stimulus characteristics
    if(compute_and_output_display_statistics) {
        calculate_display_statistics();
        //initial_display_properties.output();
        }
    assert(search_objects.size() == n_objects);
    return search_objects;
}


// CSF: vertical bars 1 X 3.5; target: red R vertical V bar; distractors: green G vertical V bar
// object size is avg (1+3.5)/2 = 2.25
void Display_creator::generate_CSF_stimulus()
{
    // fill the container of search objects, using the first n_objects locations
    // fill container with distractors
    for(int i = 0; i < n_objects; i++) {
        search_objects.push_back(Search_object("DistGrnV", jitter(all_locations[i], CSF_x_jitter, CSF_y_jitter),
            GU::Size(1., 3.5), Green_c, Rectangle_c, Vertical_c));
        }
    // if the trial is positive polarity, replace the [0] distractor with the appropriate target
    if(polarity == Polarity_e::POSITIVE) {
        // create the target, keeping a separate copy for convenience
        target = Search_object("TargetRedV", jitter(all_locations[0], CSF_x_jitter, CSF_y_jitter),
            GU::Size(1., 3.5),Red_c, Rectangle_c, Vertical_c, true);
        search_objects[0] = target;
        }
}

// COC: vertical 1 X 3.5 or horizontal 3.5 X 1 bars; target: red vertical; distractors: red horizontal, green vertical.
// limitation of distractors to two instead of three possibilities is explicit in paper - P. 1306.
// Paper is not explicit about distribution of the two kinds of distractors, especially at n = 3.
// Assuming here that there should be at least one of each kind of distractor.
// object size is avg (1+3.5)/2 = 2.25
void Display_creator::generate_COC_stimulus()
{
    // fill the container of search objects, using the first n_objects locations
    // first fill with distractors, an equal number of each type,
    // then if a positive trial, the target object will be put into the first cell of the container.
    // so first fill container with distractors; since there are two distractor types, flip a coin to
    // decide which distractor fills the first half of the container, and which fills the second half.
    // This way, the target replaces an equal number of distractors of each type over trials.
    // n_objects == 3 is special case because it is odd - fill the first two cells with a distractor,
    // the third cell with the other distractor, ensuring that at least one distractor of each type is present.

    // create skeletons for each kind of distractor
    Search_object dist1("DistRedH", GU::Size(3.5, 1.), Red_c, Rectangle_c, Horizontal_c);
    double x_jitter1 = COCH_x_jitter;
    double y_jitter1 = COCH_y_jitter;
    Search_object dist2("DistGrnV", GU::Size(1., 3.5), Green_c, Rectangle_c, Vertical_c);
    double x_jitter2 = COCV_x_jitter;
    double y_jitter2 = COCV_y_jitter;
    // now flip a coin and swap the two skeletons to decide who goes first into the container
    if(biased_coin_flip(0.5)) {
        swap(dist1, dist2);
        swap(x_jitter1, x_jitter2);
        swap(y_jitter1, y_jitter2);
        }
    // n_objects is 3, 6, 12, or 18; the "extra" distractor has been chosen at random
    // for n=3, ensure that at least one distractor of each type is chosen by using two of the first kind, one of the other.
    int n_half_objects = (n_objects == 3) ? 2 : n_objects / 2;
    for(int i = 0; i < n_half_objects; i++) {
        dist1.set_location(jitter(all_locations[i], x_jitter1, y_jitter1));
        dist1.set_ID(); // create unique object ID
        search_objects.push_back(dist1); // copy it into the container.
        }
    for(int i = n_half_objects; i < n_objects; i++) {
        dist2.set_location(jitter(all_locations[i], x_jitter2, y_jitter2));
        dist2.set_ID(); // create unique object ID
        search_objects.push_back(dist2); // copy it into the container.
        }
/*
    for(int i = 0; i < n_half_objects; i++) {
        search_objects.push_back(Search_object("DistRedH", jitter(all_locations[i], COCH_x_jitter, COCH_y_jitter),
            GU::Size(3.5, 1.), Red_c, Rectangle_c, Horizontal_c));
        }
    for(int i = n_half_objects; i < n_objects; i++)
        search_objects.push_back(Search_object("DistGrnV", jitter(all_locations[i], COCV_x_jitter, COCV_y_jitter),
            GU::Size(1., 3.5), Green_c, Rectangle_c, Vertical_c));
*/
    // if the trial is positive polarity, replace the [0] distractor with the appropriate target
    if(polarity == Polarity_e::POSITIVE) {
        // create the target, keeping a separate copy for convenience
        target = Search_object("TargRedV", jitter(all_locations[0], COCV_x_jitter, COCV_y_jitter),
            GU::Size(1., 3.5), Red_c, Rectangle_c, Vertical_c, true);
        search_objects[0] = target;
        }
}

// SHP: "digital" 2s and 5s 1.5 X 2.7; target: 2, distractors 5
// object size is avg (1.5+2.7)/2 = 2.1

void Display_creator::generate_SHP_stimulus()
{
    // fill the container of search objects, using the first n_objects locations
    // fill container with distractors
    for(int i = 0; i < n_objects; i++) {
//        search_objects.push_back(Search_object("Dist5", jitter(all_locations[i], SHP_x_jitter, SHP_y_jitter),
//            GU::Size(1.5, 2.7), Black_c, Circle_c, Nil_c));
        search_objects.push_back(Search_object("Dist5", jitter(all_locations[i], SHP_x_jitter, SHP_y_jitter),
            GU::Size(1.5, 2.7), Black_c, Digital_5_c, Vertical_c));
        }
    // if the trial is positive polarity, replace the [0] distractor with the appropriate target
    if(polarity == Polarity_e::POSITIVE) {
        // create the target, keeping a separate copy for convenience
//        target = Search_object("Targ2", jitter(all_locations[0], SHP_x_jitter, SHP_y_jitter),
//            GU::Size(1.5, 2.7), Black_c, Square_c, Nil_c);
        target = Search_object("Targ2", jitter(all_locations[0], SHP_x_jitter, SHP_y_jitter),
            GU::Size(1.5, 2.7), Black_c, Digital_2_c, Vertical_c, true);
        search_objects[0] = target;
        }
}

void Display_properties_accumulator::output() const
{
    Ostream_format_saver ofs(cout);
    // Set output to show two decimal places for floating point numbers
//    cout << fixed << setprecision(2) << endl;
    cout.setf(ios::fixed, ios::floatfield);
    cout.precision(2);
    cout << eccentricities.get_mean() << '\t'
        << spacings.get_mean() << '\t'
        << numbers_of_crowding_flankers.get_mean() << '\t'
        << dbn_of_number_of_crowding_flankers.get_max() << '\t'
        << spacings_of_crowding_flankers.get_mean() << ":\t";
    auto dbn = dbn_of_number_of_crowding_flankers.get_distribution();
    for(auto x : dbn)
        cout << '\t' << x;
    cout << endl;
}

void Display_creator::calculate_display_metrics()
{
    // calculate the display metrics for every possible object i and other object j (!= i)
    for(int i = 0; i < search_objects.size(); i++) { // i is current object
        // find and save nearest neighbor distance
        double nearest_neighbor_distance = 99.;
        for(int j = 0; j < search_objects.size(); j++) { // j is possible other object
            if(j != i) { // j is other object
                double eccentricity = GU::cartesian_distance(search_objects[i].location, search_objects[j].location);
                if(eccentricity < nearest_neighbor_distance)
                    nearest_neighbor_distance = eccentricity;
                } // end of possible other object j
            } // end of possible other object loop
        search_objects[i].nearest_neighbor_distance = nearest_neighbor_distance;
        } // end of current object i
}

void Display_creator::reset_display_statistics()
{
    initial_display_properties.reset();
    whole_display_properties.reset(); // why reset here?
}

// let each object be the current object, calculate the eccentricity to each other object (target),
// and for that target object calculate the number of flanking objects within the critical spacing (Bouma)
// calculate similar statistics for eccentricity from fixation point
void Display_creator::calculate_display_statistics()
{
    initial_display_properties.n_trials++;
	initial_display_properties.number_of_objects.update(search_objects.size());
    whole_display_properties.n_trials++;
	whole_display_properties.number_of_objects.update(search_objects.size());

    // calculating the initial display statistics is a special case because the fixation point is not a search_object
    GU::Point current_fixation_location(0., 0.);
    for(int j = 0; j < search_objects.size(); j++) {
        // each object j is possible target for next fixation
        calculate_possible_next_fixation_statistics(current_fixation_location, j, initial_display_properties);
        }
    
    // calculate the display statistics for every possible current fixation i and next fixation j (!= i) and possible flanker k (!= j)
    for(int i = 0; i < search_objects.size(); i++) { // object i is current point of fixation
        GU::Point current_fixation_location = search_objects[i].location;
        for(int j = 0; j < search_objects.size(); j++) {
            if(j != i) { // object j is a possible target
                calculate_possible_next_fixation_statistics(current_fixation_location, j, whole_display_properties);
                } // end of possible target j
            } // end of possible target loop
        } // end of current fixation point i
}

void Display_creator::calculate_possible_next_fixation_statistics(GU::Point current_assumed_fixation_point, int possible_next_fixation_index, Display_properties_accumulator& properties_accumulator)
{
    GU::Point possible_next_fixation_location = search_objects[possible_next_fixation_index].location;
    double eccentricity = GU::cartesian_distance(current_assumed_fixation_point, possible_next_fixation_location);
    properties_accumulator.eccentricities.update(eccentricity);
    properties_accumulator.dbn_of_eccentricities.update(eccentricity);
    int n_crowders = 0; // how many flanker objects are crowding this possible target possible_next_fixation_index
    for(int k = 0; k < search_objects.size(); k++) {
        if(k != possible_next_fixation_index) { // k is a possible flanker of target possible_next_fixation_index (current point of fixation could be a flanker)
            double spacing = GU::cartesian_distance(possible_next_fixation_location, search_objects[k].location);
            properties_accumulator.spacings.update(spacing);
            if(spacing <= bouma_fraction_c * eccentricity) {
                // flanker k is crowding target j
                n_crowders++;
                // if this flanker is within crowding distance, update the average for that spacing
                properties_accumulator.spacings_of_crowding_flankers.update(spacing);
                } // end of crowding flanker
            } // end of possible flanker k
        }
    properties_accumulator.numbers_of_crowding_flankers.update(n_crowders);
    properties_accumulator.dbn_of_number_of_crowding_flankers.update(n_crowders);
    properties_accumulator.mean_n_crowders_for_eccentricity.update(eccentricity, n_crowders);
}


void Display_creator::output_display_statistics() const
{
    output_display_properties_accumulator("Initial display properties", initial_display_properties);
    output_display_properties_accumulator("Whole display properties", whole_display_properties);
}

void Display_creator::output_display_properties_accumulator(const char* message, const Display_properties_accumulator& acc) const
{
    cout << message << endl;
//	cout << "Ntrials" << '\t' << "NObjs" << '\t' <<
	cout << "Ntrials" << '\t' <<
		"NObMin" << '\t' << "NObM" << '\t' << "NObMx" << '\t' <<
		"EccMin" << '\t' << "EccM" << '\t' << "EccMx" << '\t' <<
		"SpaMin" << '\t' << "SpaM" << '\t' << "SpaMx" << '\t' <<
		"NCFMin" << '\t' << "NCFM" << '\t' << "NCFMx" << '\t' <<
		"SCFMin" << '\t' << "SCFM" << '\t' << "SCFMx" << endl;
	cout << acc.n_trials << '\t' <<
		acc.number_of_objects.get_min() << '\t' << acc.number_of_objects.get_mean() << '\t' << acc.number_of_objects.get_max() << '\t' <<
		acc.eccentricities.get_min() << '\t' << acc.eccentricities.get_mean() << '\t' << acc.eccentricities.get_max() << '\t' <<
		acc.spacings.get_min() << '\t' << acc.spacings.get_mean() << '\t' << acc.spacings.get_max() << '\t' <<
		acc.numbers_of_crowding_flankers.get_min() << '\t' << acc.numbers_of_crowding_flankers.get_mean() << '\t' << acc.numbers_of_crowding_flankers.get_max() << '\t' <<
		acc.spacings_of_crowding_flankers.get_min() << '\t' << acc.spacings_of_crowding_flankers.get_mean() << '\t' << acc.spacings_of_crowding_flankers.get_max() << endl;
	// output distributions
	cout << "EccDbn"  << '\t' << acc.dbn_of_eccentricities.get_min() << '\t' << acc.dbn_of_eccentricities.get_max() << '\t' << acc.dbn_of_eccentricities.get_n();
	for(int i = 0; i < acc.dbn_of_eccentricities.get_n_bins(); i++)
		cout << '\t' << acc.dbn_of_eccentricities[i];
	cout << endl;
	cout << "NCFDbn"  << '\t' << acc.dbn_of_number_of_crowding_flankers.get_min() << '\t' << acc.dbn_of_number_of_crowding_flankers.get_max() << '\t' << acc.dbn_of_number_of_crowding_flankers.get_n();
	for(int i = 0; i < acc.dbn_of_number_of_crowding_flankers.get_n_bins(); i++)
		cout << '\t' << acc.dbn_of_number_of_crowding_flankers[i];
	cout << endl;
	cout << "MNCEcc"  << '\t' << acc.mean_n_crowders_for_eccentricity.get_min() << '\t' << acc.mean_n_crowders_for_eccentricity.get_max() << '\t' << acc.mean_n_crowders_for_eccentricity.get_n();
	for(int i = 0; i < acc.mean_n_crowders_for_eccentricity.get_n_bins(); i++)
		cout << '\t' << acc.mean_n_crowders_for_eccentricity[i].get_mean();
	cout << endl;
}


GU::Point jitter(GU::Point point, double x_jitter, double y_jitter)
{
    return GU::Point(uniform_random_variable(point.x, x_jitter), uniform_random_variable(point.y, y_jitter));
}




------------------------------------------------------------
// Geometry.h
#ifndef GEOMETRY_H
#define GEOMETRY_H

#define USE_ASSERT_THROW

#ifdef USE_ASSERT_THROW
#include "Assert_throw.h"
#define ASSERT Assert
#endif

#ifdef USE_CASSERT
#include <cassert>
#define ASSERT assert
#endif

#include "Point.h"
#include <iosfwd>
#include <vector>
#define _USE_MATH_DEFINES  // from http://stackoverflow.com/questions/1727881/how-to-use-the-pi-constant-in-c
#include <cmath>

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#ifndef M_PIl
#define M_PIl (3.14159265358979323846264338327950288)
#endif


// calculate a value for pi - this has internal linkage
//const double GU_pi = 2. * std::atan2(1., 0.);
const double GU_pi = M_PI;

namespace Geometry_Utilities {
using std::operator<<;	// bring in global declarations
struct Point;
struct Size;
struct Cartesian_vector;
struct Polar_vector;

// Output operators
std::ostream& operator<< (std::ostream& os, const Point& p);
std::ostream& operator<< (std::ostream& os, const Size& s);
std::ostream& operator<< (std::ostream& os, const Cartesian_vector& cv);
std::ostream& operator<< (std::ostream& os, const Polar_vector& pv);

/*
This set of simple classes is used to compute positions and directions in the plane, 
in both cartesian and polar coordinates, and also with respect to line segments
and polygon regions.

Most classes are defined with "struct" to make all members public by default.

These classes make no assumptions about units of measurement of distance. Angles
can be specified in radians, but radians can be converted to and from 
trigonometry degrees in which 0 degrees corresponds to (x > 0, y = 0), 
and 90 degrees corresponds to (x = 0, y > 0).

Point is a set of (x, y) coordinates. 

A Cartesian_vector is (delta_x, delta_y) - a displacement in Cartesian coordinates.

A Polar_vector is (r, theta) - a displacement in polar coordinates using radians.

Various overloaded operators support computations of positions and directions.

A Line_segment represents a line that passes through two Points in 
both parametric and general form. Functions are provided (mostly inline for speed)
for computing distances of a Point from the line, and intersections between lines.
There are two definitions supported for distance of a point from the line:
	1. The distance from the line extended infinitely past the endpoints of the segment.
	2. The distance from the line segment only, which is either:
		a. the length of the line that both passes through the point and is perpendicular
		to the line segment and intersects the line segment between the endpoints. 
		b. the distance from the closest endpoint if the intersecting perpendicular does
		not intersect between the endpoints.

The Polygon class represents a polygon as a series of Points for vertices that make up
the endpoints of a set of line segments. The last Point wraps around to the first point
to define the last line segment making up the polygon. A distance_inside function
is used to compute whether a Point is inside the polygon and the distance of the Point
from the nearest line segment.

*/

// calculate a value for pi - this should have internal linkage
// const double pi = 2. * std::atan2(1., 0.);
//extern const double pi;

// angle units conversion functions
double to_radians(double theta_d);
double to_degrees(double theta_r);

// use these functions for things like visual angle per pixel
double degrees_subtended(double size_measure, double distance_measure);
double degrees_subtended_per_unit(double units_per_measure, double distance_measure);
double units_per_degree_subtended(double units_per_measure, double distance_measure);

// forward class declarations
struct Cartesian_vector;
struct Polar_vector;


/* 
// A Point contains an (x, y) pair to represent coordinates
struct Point
{
	double x;
	double y;
	
	Point (double in_x = 0., double in_y = 0.) :
		x(in_x), y(in_y)
		{}

	// compare two Points
	bool operator== (const Point& rhs) const
		{return (x == rhs.x && y == rhs.y);}
 
	bool operator!= (const Point& rhs) const
		{return !(*this == rhs);}
		
	bool operator< (const Point& rhs) const
		{return (x == rhs.x) ? (y < rhs.y) : (x < rhs.x);}
		
	bool operator<= (const Point& rhs) const
		{return (*this < rhs) || (*this == rhs);}
		
	bool operator> (const Point& rhs) const
		{return !(*this <= rhs);}

	bool operator>= (const Point& rhs) const
		{return !(*this < rhs);}
};
*/

/* Size */
// A Size contains an (h, v) pair to represent the horizontal
// and vertical size of an object, which is different from
// a Point or a Cartesian_vector.
struct Size
{
	double h;
	double v;
	
	explicit Size (double in_h = 0., double in_v = 0.) :
		h(in_h), v(in_v)
		{}

	// compare two Sizes
	bool operator== (const Size& rhs) const
		{
			return (h == rhs.h && v == rhs.v);
		}
 
	bool operator!= (const Size& rhs) const
		{
			return (h != rhs.h || v != rhs.v);
		}
};


// return the distance between two Points
double cartesian_distance (const Point& p1, const Point& p2);

/* Cartesian_vector */
// A Cartesian_vector contains an x, y displacement
struct Cartesian_vector
{
	double delta_x;
	double delta_y;

	Cartesian_vector (double in_delta_x = 0., double in_delta_y = 0.) : 
		delta_x(in_delta_x), delta_y(in_delta_y)
	{}

	// construct a Cartesian_vector from two Points,
	// showing the vector from p1 to p2 
	// that is, p1 + cv => p2
	Cartesian_vector(const Point& p1, const Point& p2);

	// construct a Cartesian_vector from a Polar_vector
	Cartesian_vector(const Polar_vector& pv);
};


/* Polar_vector */
// Polar_vector describes a displacement in terms of polar coordinates
// with angle in radians
struct Polar_vector
{
	double r;
	double theta;

	Polar_vector (double in_r = 0., double in_theta = 0.) : 
		r(in_r), theta(in_theta)
	{}

	// construct a Polar_vector from two Points,
	// showing the vector from p1 to p2 
	// that is, p1 + pv => p2
	Polar_vector(const Point& p1, const Point& p2);

	// construct a Polar_vector from a Cartesian_vector
	Polar_vector(const Cartesian_vector& cv);

};




// *** Overloaded Operators ***
// For ease in documentation, all overloaded operators are non-member functions

// Subtract two Points to get a Cartesian_vector
// p2's components are subtracted from p1
Cartesian_vector operator- (const Point& p1, const Point& p2);

// Add a Point and a Cartesian_vector to get the displaced Point
Point operator+ (const Point& p, const Cartesian_vector& cv);
Point operator+ (const Cartesian_vector& cv, const Point& p);
	
// Add a Point and a Polar_vector to get the displaced Point
Point operator+ (const Point& p, const Polar_vector& pv);
Point operator+ (const Polar_vector& pv, const Point& p);

// Adding or subtracting two Cartesian_vectors adds or subtracts the components
Cartesian_vector operator+ (const Cartesian_vector& cv1, const Cartesian_vector& cv2);
Cartesian_vector operator- (const Cartesian_vector& cv1, const Cartesian_vector& cv2);

// divide a Cartesian_vector by a double: divide each component by the double
Cartesian_vector operator/ (const Cartesian_vector& cv, double d);
Cartesian_vector operator/ (double d, const Cartesian_vector& cv);

// divide a Polar_vector by a double: divide r component by the double
Polar_vector operator/ (const Polar_vector& pv, double d);
Polar_vector operator/ (double d, const Polar_vector& pv);

// multiply a Cartesian_vector by a double: multiply each component by the double
Cartesian_vector operator* (const Cartesian_vector& cv, double d);
Cartesian_vector operator* (double d, const Cartesian_vector& cv);

// multiply a Polar_vector by a double: multiply r component by the double
Polar_vector operator* (const Polar_vector& pv, double d);
Polar_vector operator* (double d, const Polar_vector& pv);

/*
The Line_segment class stores two endpoints and the terms of the general form,
and precomputes for speed some values used in various calculations.

In the parametric form, line is specified in terms of its endpoints, p1 and p2, 
the differences dx = p2.x - p1.x, and dy = p2.y - p1.y, and a parameter t
that describes where on the line a particular point is, starting from p1.

A point is on the line if its (x, y) coordinates can be specified as
    x = p1.x + t * dx
    y = p1.y + t * dy 
    for some value of t.
Note that p1 corresponds to t = 0, p2 to t = 1. So for a given (x, y), a value
of t between 0 and 1 means that the point lies between the endpoints. A value
of t < 0 means that t comes "before" p1


In the general form, a point (x, y) is on the line if it satisfies the equation
	A * x + B * y + C == 0, where
	
	len = distance between p1 and p2
	A = -dy/len;
	B = dx/len;
	C = -(-dy*x1 + dx*y1)/len

For speed, most functions are inline, and dx, dy, len, and c (numerator of C) 
are computed and saved when Line_segment is initialized. 

Note that x and y values are doubles, and if they have fractional values, some
functions and tests that assume exact values will not be accurate. 
But if all coordinates have integral values, then functions such as is_horizontal 
should be correct.
*/

class Line_segment {
public:
	Line_segment() : // default is all values zero - normally not used.
		dx(0.), dy(0.), c(0.), len(0.)
		{}
	Line_segment(Point in_p1, Point in_p2);
	
	Point get_p1() const
		{return p1;}

	Point get_p2() const
		{return p2;}

	// return the center of the bounding box
	Point get_center() const
		{return center;}
	// return the Size of the bounding box
	Size get_size() const
		{return size;}

	double get_dx() const
		{return dx;}

	double get_dy() const
		{return dy;}

	double get_A() const
		{return -dy/len;}

	double get_B() const
		{return dx/len;}

	double get_C() const
		{return c/len;}

	double get_length() const
		{return len;}

	bool is_horizontal() const
		{return dy == 0.;}

	bool is_vertical() const
		{return dx == 0.;}
	
	// These functions treat the segment as an infinite line
	// Return true if p is on the infinite line described by the segment
	bool is_on_infinite_line(Point p) const
		{
			if((-dy * p.x + dx * p.y + c) == 0)
				return true;
			return false;
		}

	// Compute the closest point on the infinite line
	Point closest_point_on_infinite_line(Point p) const
		{return point_on_line(parameter(p));}

	// Return the distance from the line to a point
	double distance_from_infinite_line(Point p) const;

	// Return the parameter value for the point. 
	// The value will be between 0 and 1 if the point is between the endpoints.
	// Otherwise, t corresponds to the position of the point along the line.
	// t < 0 means p is "before" p1, t > 1, "after" p2.
	double parameter(Point p) const
		{return (-(p1.x - p.x) * dx - (p1.y - p.y) * dy) / (len * len);}

	// Compute t given a value for x
	double parameter_given_x(double x) const
		{
			ASSERT (dx != 0.); // vertical line is indeterminate!
			return (x - p1.x) / dx;
		}

	// Return x of point on the line specified by t
	double x_given_parameter(double t) const
		{return p1.x + t * dx;}

	// Compute t given a value for y
	double parameter_given_y(double y) const
		{
			ASSERT(dy != 0.);	// horizontal line is indeterminate!
			return (y - p1.y) / dy;
		}

	// Return y of point on the line specified by t
	double y_given_parameter(double t) const
		{return p1.y + t * dy;}

	// Return the point on the line specified by t
	Point point_on_line(double t) const
		{return Point(p1.x + t * dx, p1.y + t * dy);}

	// Return the closest point on the line segment to p;
	// This is either a point on the line, or the closest endpoint
	Point closest_point_on_segment(Point p) const
		{
			double t = parameter(p);
			if (t <= 0.)
				return p1;
			else if (t >= 1.)
				return p2;
			else
				return point_on_line(t);
		}

	// Return the distance from p to the closest point on the line segment
	double distance_from_segment(Point p) const
		// {return distance(p, closest_point_on_segment(p));}
		{
			double t = parameter(p);
			if (t <= 0.)
				return cartesian_distance(p, p1);
			else if (t >= 1.)
				return cartesian_distance(p, p2);
			else {
				return distance_from_infinite_line(p);
				}
		}

private:
	Point p1;	// the "start" point
	Point p2;	// the "stop" point
	double dx;	// (x2 - x1) is -a, A is -a/len
	double dy;	// (y2 - y1) is b, B is b/len
	double c;	// C is c/len
	double len;	// length of the line segment

	Point center;	// center of the bounding box
	Size size;	// size of the bounding box

};

std::ostream& operator<< (std::ostream& os, const Line_segment& ls);

// Given a point and a rectangle, return true if the point is inside the rectangle,
// false if not.
bool is_point_inside_rectangle(Point p, Point rect_loc, Size rect_size);

// Given a line and a rectangle, compute and return the line segment that is the line clipped to the rectangle.
Line_segment clip_line_to_rectangle(const Line_segment& line, Point rect_loc, Size rect_size);

// Given a line from a point through the center of a rectangle, calculate
// the line that goes from the closest point of intersection on the rectangle
// to the center of the rectangle. Use for e.g. Fitts ID calculations.
bool compute_center_intersecting_line(const Line_segment& start_to_center, Size rect_size, Line_segment& clipped_line);

// Compute the closest distance from p to the rectangle given by center, size
double closest_distance(Point p, Point rect_center, Size rect_size);

// The polygon class stores a vector of points corresponding to adjacent vertices
// of the polygon - the last point "wraps" to the first point to close the polygon.
class Polygon {
public:
	// an empty polygon has default size and center
	Polygon ()	: cache_good(true) 
		{}
	// Can initialize with a vector of points if desired
	Polygon(const std::vector<Point>& in_vertices) :
		vertices(in_vertices), cache_good(false)
		{}
	// create the polygon by adding points in order
	void add_vertex(Point p)
		{vertices.push_back(p); cache_good = false;}
	// empty the polygon
	void clear()
		{vertices.clear(); cache_good = false;}
	// supply a const reference to the vertices
	const std::vector<Point>& get_vertices() const
		{return vertices;}
	// return the center of the bounding box
	Point get_center() const;
	// return the Size of the bounding box
	Size get_size() const;
	// compute the distance of p from the edge of the polygon. 
	// A negative value means that p is outside the polygon, positive means p is inside
	double distance_inside(const Point p) const;

private:
	std::vector<Point> vertices;
	// cached values are recomputed if points changed since last calculation
	mutable bool cache_good;	// false if size and center need to be computed anew
	mutable Point center;
	mutable Size size;
	void recompute() const;	// const since called from const functions
};

} // end namespace


#endif

------------------------------------------------------------
// Geometry.cpp
/* 
Geometry.cpp implementation file
See Geometry.h for comments
*/
#include "Geometry.h"

#include <iostream>
#include <vector>
#include <cmath>

using std::ostream;
using std::vector;
using std::size_t;
using std::sqrt;	using std::cos;	using std::sin;	using std::atan2;	using std::fabs;


namespace Geometry_Utilities {



// *** Member function definitions ***

// return the distance between two Points
double cartesian_distance (const Point& p1, const Point& p2)
{
	double xd = p2.x - p1.x;
	double yd = p2.y - p1.y;
	double d = sqrt(xd * xd + yd * yd);
	return d;
}

// Cartesian_vector members
// construct a Cartesian_vector from two Points,
// showing the vector from p1 to p2 
// that is, p1 + cv => p2
Cartesian_vector::Cartesian_vector(const Point& p1, const Point& p2)
{
	delta_x = p2.x - p1.x;
	delta_y = p2.y - p1.y;
}

// construct a Cartesian_vector from a Polar_vector
Cartesian_vector::Cartesian_vector(const Polar_vector& pv)
{
	delta_x = pv.r * cos(pv.theta);
	delta_y = pv.r * sin(pv.theta);
}

// Polar_vector members
// construct a Polar_vector from a Cartesian_vector
// Theta ranges from [0 to 2*pi)
/*
pi = 3.14159 2*pi = 6.28319
Polar_vector(Point(0,0), Point(+1, 0.))	P<1, 0>
Polar_vector(Point(0,0), Point(+1, 0.01))	P<1.00005, 0.00999967>
Polar_vector(Point(0,0), Point(+1, +1))	P<1.41421, 0.785398>
Polar_vector(Point(0,0), Point(0., +1))	P<1, 1.5708>
Polar_vector(Point(0,0), Point(-1, +1))	P<1.41421, 2.35619>
Polar_vector(Point(0,0), Point(-1, 0.))	P<1, 3.14159>
Polar_vector(Point(0,0), Point(-1, -1))	P<1.41421, 3.92699>
Polar_vector(Point(0,0), Point(0., -1))	P<1, 4.71239>
Polar_vector(Point(0,0), Point(+1, -1))	P<1.41421, 5.49779>
Polar_vector(Point(0,0), Point(+1, -0.5))	P<1.11803, 5.81954>
Polar_vector(Point(0,0), Point(+1, -0.1))	P<1.00499, 6.18352>
Polar_vector(Point(0,0), Point(+1, -0.01))	P<1.00005, 6.27319>
Polar_vector(Point(0,0), Point(+1, -0.001))	P<1, 6.28219>
*/
Polar_vector::Polar_vector(const Cartesian_vector& cv)
{
	r = sqrt ((cv.delta_x * cv.delta_x) + (cv.delta_y * cv.delta_y));
	// atan2 will return neg angle for Quadrant III, IV.
	theta = atan2 (cv.delta_y, cv.delta_x);
	if (theta < 0.)
		theta = 2. * GU_pi + theta; // normalize theta positive
}

// construct a Polar_vector from two Points,
// showing the vector from p1 to p2 
// that is, p1 + pv => p2
Polar_vector::Polar_vector(const Point& p1, const Point& p2)
{
	Polar_vector pv (Cartesian_vector(p1, p2));
	r = pv.r;
	theta = pv.theta;
}




// *** Overloaded Operators ***
// For ease in documentation, all overloaded operators are non-member functions

// Subtract two Points to get a Cartesian_vector
// p2's components are subtracted from p1
Cartesian_vector operator- (const Point& p1, const Point& p2)
{
	return Cartesian_vector(p1.x - p2.x, p1.y - p2.y);
}

// Add a Point and a Cartesian_vector to get the displaced Point
Point operator+ (const Point& p, const Cartesian_vector& cv)
{
	return Point(p.x + cv.delta_x, p.y + cv.delta_y);
}

Point operator+ (const Cartesian_vector& cv, const Point& p)
{
	return p + cv;
}
	
// Add a Point and a Polar_vector to get the displaced Point
Point operator+ (const Point& p, const Polar_vector& pv)
{
	Cartesian_vector cv (pv);
	return cv + p;
}

Point operator+ (const Polar_vector& pv, const Point& p)
{
	return p + pv;
}

// Adding or subtracting two Cartesian_vectors adds or subtracts the components
Cartesian_vector operator+ (const Cartesian_vector& cv1, const Cartesian_vector& cv2)
{
	return Cartesian_vector(cv1.delta_x + cv2.delta_x, cv1.delta_y + cv2.delta_y);
}

Cartesian_vector operator- (const Cartesian_vector& cv1, const Cartesian_vector& cv2)
{
	return Cartesian_vector(cv1.delta_x - cv2.delta_x, cv1.delta_y - cv2.delta_y);
}

// divide a Cartesian_vector by a double: divide each component by the double
Cartesian_vector operator/ (const Cartesian_vector& cv, double d)
{
	return Cartesian_vector(cv.delta_x / d, cv.delta_y / d);
}

Cartesian_vector operator/ (double d, const Cartesian_vector& cv)
{
	return cv / d;
}

// divide a Polar_vector by a double: divide r component by the double
Polar_vector operator/ (const Polar_vector& pv, double d)
{
	return Polar_vector(pv.r / d, pv.theta);
}
Polar_vector operator/ (double d, const Polar_vector& pv)
{
	return pv / d;
}

// multiply a Cartesian_vector by a double: multiply each component by the double
Cartesian_vector operator* (const Cartesian_vector& cv, double d)
{
	return Cartesian_vector(cv.delta_x * d, cv.delta_y * d);
}

Cartesian_vector operator* (double d, const Cartesian_vector& cv)
{
	return cv * d;
}

// multiply a Polar_vector by a double: multiply r component by the double
Polar_vector operator* (const Polar_vector& pv, double d)
{
	return Polar_vector(pv.r * d, pv.theta);
}

Polar_vector operator* (double d, const Polar_vector& pv)
{
	return pv * d;
}

// Output operators
// output a Point as "(x, y)"
ostream& operator<< (ostream& os, const Point& p)
{
	os << '(' << p.x << ", " << p.y << ')';
	return os;
}
	
// output a Size as "[h, v]"
ostream& operator<< (ostream& os, const Size& s)
{
	os << '[' << s.h << ", " << s.v << ']';
	return os;
}

// output a Cartesian_vector as "<x, y>"
ostream& operator<< (ostream& os, const Cartesian_vector& cv)
{
	os << '<' << cv.delta_x << ", " << cv.delta_y << '>';
	return os;
}

// output a Polar_vector as "P<r, theta>"
ostream& operator<< (ostream& os, const Polar_vector& pv)
{
	os << "P<" << pv.r << ", " << pv.theta << '>';
	return os;
}

/***** Utility function definitions *****/
// There are 2pi radians in 360 degrees
double to_radians (double theta_d)
{
	return 2. * GU_pi  * (theta_d / 360);
}

double to_degrees (double theta_r)
{
	const double conversion_factor = 180./GU_pi;
	double temp = theta_r * conversion_factor;
	return temp;
}

// size_measure and distance_measure must be in the same units (e.g. pixels, inches)
// if signed, atan2 uses these to determine the quadrants involved.
double degrees_subtended(double size_measure, double distance_measure)
{
	double radians = 2. * atan2(size_measure, 2. * distance_measure);
	return to_degrees(radians);
}

// I am not sure of the purpose of this definition - seems roundabout
double degrees_subtended_per_unit(double units_per_measure, double distance_measure)
{
    return degrees_subtended(1.0/units_per_measure, distance_measure);
}

double units_per_degree_subtended(double units_per_measure, double distance_measure)
	{return 1./degrees_subtended_per_unit(units_per_measure, distance_measure);}

/*
The Line_segment class stores two endpoints and the terms of the general form,
and precomputes for speed some values used in various calculations.
*/
Line_segment::Line_segment(Point in_p1, Point in_p2) :
	p1(in_p1), p2(in_p2)
{
	dx = p2.x - p1.x;
	dy = p2.y - p1.y;
	len = sqrt(dx * dx + dy * dy);
	c = -(-dy * p1.x + dx * p1.y);
	
	center.x = (p1.x + p2.x) / 2.;
	center.y = (p1.y + p2.y) / 2.;

	size.h = fabs(dx);
	size.v = fabs(dy);

//	center.x = p1. x + (dx > 0.) ? size.h / 2.0 : - size.h / 2.0;
//	center.y = p1. y + (dy > 0.) ? size.v / 2.0 : - size.v / 2.0;
}

// Return the distance from the line to a point
double Line_segment::distance_from_infinite_line(Point p) const
{
	return fabs(-dy * p.x + dx * p.y + c) / len;
}

// output a Line_segment as "p1 -> p2"
ostream& operator<< (ostream& os, const Line_segment& ls)
{
	os << ls.get_p1() << " -> " << ls.get_p2();
	return os;
}

// Given a point and a rectangle specified as the center of the rectangle, and its size
// return true if the point is inside the rectangle, false if not.
bool is_point_inside_rectangle(Point p, Point rect_loc, Size rect_size)
{
    double xmin = rect_loc.x - rect_size.h / 2.;
    double xmax = rect_loc.x + rect_size.h / 2.;
    if(p.x < xmin || p.x > xmax)
        return false;
    double ymin = rect_loc.y - rect_size.v / 2.;
    double ymax = rect_loc.y + rect_size.v / 2.;
    if(p.y < ymin || p.y > ymax)
        return false;
    return true;
}



// Given a line and a rectangle, compute and return the line segment that is the line clipped to the rectangle.
// Based on the Cohen-Sutherland clipping algorithm (see Foley, van Damm, Feiner, & Hughes).

Line_segment clip_line_to_rectangle(const Line_segment& line, Point rect_loc, Size rect_size)
{
	double xmin = rect_loc.x;
	double xmax = rect_loc.x + rect_size.h;
	double ymin = rect_loc.y;
	double ymax = rect_loc.y + rect_size.v;
	Point p1 = line.get_p1();
	Point p2 = line.get_p2();
	double dx = line.get_dx();
	double dy = line.get_dy();
	
	double txmax = line.parameter_given_x(xmax);
	double txmin = line.parameter_given_x(xmin);
	double tymax = line.parameter_given_y(ymax);
	double tymin = line.parameter_given_y(ymin);
	
	Point cp1, cp2;
	// line goes from top-right (p1) towards bottom left (p2) (updside-down y)	
	if(dx < 0. && dy > 0.) {
		// choose cp1
		if(txmax < 0. && tymin < 0.) {
			// xmax and ymin are before p1, so use p1
			cp1 = p1;
			}
		else if(txmax < 0. && tymin > 0. && tymin < 1.0) {
			// xmax is before p1, clip from tymin (upside-down, remember!)
			cp1 = line.point_on_line(tymin);
			}
		else if(tymin < 0. && txmax > 0. && txmax < 1.0) {
			// ymin is before p1, clip from txmax
			cp1 = line.point_on_line(txmax);
			}
		// choose cp2
		if(txmin > 1.0 && tymax > 1.0) {
			// xmin and ymax are "after" p2, use p2
			cp2 = p2;
			}
		else if(tymax > 1.0 && txmin > 0. && txmin < 1.0) {
			// ymax is "after" p2, clip to txmin
			cp2 = line.point_on_line(txmin);
			}
		else if(txmin > 1.0 && tymax > 0. && tymax < 1.0) {
			// xmin is "after" p2, clip to tymax
			cp2 = line.point_on_line(tymax);
			}
		}
	// line goes from bottom-left towards top-right		
	else if(dx > 0. && dy < 0.) {
		// choose cp1
		if(txmin < 0. && tymax < 0.) {
			// xmin and ymax are before p1, so use p1
			cp1 = p1;
			}
		// choose cp2
		if(txmax > 1. && tymin > 0. && tymin < 1.0) {
			// xmax is after p2, but not ymin, so use tymin
			cp2 = line.point_on_line(tymin);
			}
		else if(txmax > 1.0 && tymin > 1.0) {
			// both xmax and ymin are after p2, so use p2
			cp2 = p2;
			}
		}
		
	// line goes from top-left(p1) towards bottom-right (p2)
	else if(dx > 0. && dy > 0.) {
		// choose cp1
		if(txmin < 0. && tymin < 0.) {
			// xmin and ymin are before p1, so use p1
			cp1 = p1;
			}
		// choose cp2
		if(tymax > 1.0 && txmax > 1.0) {
			// ymax and xmax are both after p2, so use p2
			cp2 = p2;
			}
		else if(tymax > 1. && txmax > 0. && txmax < 1.0) {
			// ymax is after p2, but not xmax, so use txmax
			cp2 = line.point_on_line(txmax);
			}
		else if(txmax > 1. && tymax > 0. && tymax < 1.0) {
			// xmax is after p2, but not ymax, so use tymax
			cp2 = line.point_on_line(tymax);
			}
		else if(txmax > 0. && txmax < 1.0 && tymax > 0. && tymax < 1.0 && txmax < tymax) {
			// both xmax and ymax are before p2, but xmax is closer to p1 than ymax
			cp2 = line.point_on_line(txmax);
			}
		else if(txmax > 0. && txmax < 1.0 && tymax > 0. && tymax < 1.0 && tymax < txmax) {
			// both xmax and ymax are before p2, but ymax is closer to p1 than xmax
			cp2 = line.point_on_line(tymax);
			}
		}

	
	Line_segment clipped_line(cp1, cp2);

	return clipped_line;

}


/*
Calculate the line that intersects a rectangle through its center, using a specialization 
of the Cohen-Sutherland clipping algorithm (see Foley, van Damm, Feiner, & Hughes).

The start-to-center line goes from the start point P0 to the end point P1, 
which is the center of the rectangular target region, whose boundaries are 
min_edges.x, max_edges.x, min_edges.y, max_edges.y.

Since P1 is always inside the rectangle, we check only for whether P0 is also inside,
and then which edge of the rectangle the line crosses, and compute the point of
intersection.

clipped_line is the line ending in the center that intersects the edge of rectangle.
return true if starting point is outside and d and s values are valid; false otherwise
*/

bool compute_center_intersecting_line(const Line_segment& start_to_center, const Size rect_size, Line_segment& clipped_line)
{
	// the start of the line is p0, the starting point, and will get moved to points of intersection during the computation
	Point p0(start_to_center.get_p1());

	// the end of the line is p1, the center of the target, and stays fixed.
	const Point p1(start_to_center.get_p2()); 

	// max and min_edges are the boundaries of the target rectangle
	// max edges are the top right coordinates
	Point max_edges(p1.x + rect_size.h/2, p1.y + rect_size.v/2);
	// min edges are the lower left coordinates
	Point min_edges(p1.x - rect_size.h/2, p1.y - rect_size.v/2);

	// Where does the line p0-p1 intersect the rectangle? Calculate intersection accordingly.
	// if the point is not outside, then start point must be inside the rectangle, 
	// and no movement should be made
	
	Point p; // point of intersection
	bool top, bottom, right, left;
	bool first = true;

	while(true) {
	top = bottom = left = right = false;
	// calculate where p0 is outside
	if (p0.y > max_edges.y) {
		top = true;
		}
	else if (p0.y < min_edges.y) {
		bottom = true;
		}
	if (p0.x > max_edges.x) {
		right = true;
		}
	else if (p0.x < min_edges.x) {
		left = true;
		}
	// if p0 is not outside, we are done
	// If it was the first value of p0, then starting point was inside, and result is invalid
	// otherwise stop and compute the final result
	// if p0 is still outside, continue the computation with the new value
	if (!(top || bottom || right || left)) {
		if(first)
			return false;
		else
			break;
		}
	else
		first = false;
	// Where does the line intersect the rectangle?
	// On the top
	if (top) {
		p.x = p1.x + (p0.x - p1.x) * (max_edges.y - p1.y)/(p0.y - p1.y);
		p.y = max_edges.y;
		}

	// On the bottom
	else if (bottom) {
		p.x = p1.x + (p0.x - p1.x) * (min_edges.y - p1.y)/(p0.y - p1.y);
		p.y = min_edges.y;
		}

	// On the right
	else if (right) {
		p.x = max_edges.x;
		p.y = p1.y + (p0.y - p1.y) * (max_edges.x - p1.x)/ (p0.x - p1.x);
		}

	// On the left
	else if (left) {
		p.x = min_edges.x;
		p.y = p1.y + (p0.y - p1.y) * (min_edges.x - p1.x)/ (p0.x - p1.x);
		}
	
	// move p0 to the point of intersection and repeat
	p0 = p;
	}	// main loop
	// return the line segment between the intersection and the center
	clipped_line = Line_segment(p0, p1);
	return true;
}

// Compute the closest distance from p to the rectangle given by center, size
double closest_distance(Point p, Point rect_center, Size rect_size)
{
	// max and min_edges are the boundaries of the target rectangle
	// tr edges are the top right coordinates
	Point tr_edges(rect_center.x + rect_size.h/2, rect_center.y + rect_size.v/2);
	// ll edges are the lower left coordinates
	Point ll_edges(rect_center.x - rect_size.h/2, rect_center.y - rect_size.v/2);
	
	double x_ll_distance = fabs(p.x - ll_edges.x);
	double x_tr_distance = fabs(p.x - tr_edges.x);
	double y_ll_distance = fabs(p.y - ll_edges.y);
	double y_tr_distance = fabs(p.y - tr_edges.y);
	
	Point vseg_p1, vseg_p2, hseg_p1, hseg_p2;
		
	if(x_ll_distance < x_tr_distance) {
		// closest vertical edge is at left
		vseg_p1.x = vseg_p2.x = ll_edges.x;
		vseg_p1.y = ll_edges.y;
		vseg_p2.y = tr_edges.y;
		}
	else {
		// closest vertical edge is at right
		vseg_p1.x = vseg_p2.x = tr_edges.x;
		vseg_p1.y = ll_edges.y;
		vseg_p2.y = tr_edges.y;
		}
	
	if(y_ll_distance < y_tr_distance) {
		// closest horizontal edge is the lower one
		hseg_p1.y = hseg_p2.y = ll_edges.y;
		hseg_p1.x = ll_edges.x;
		hseg_p2.x = tr_edges.x;
		}
	else {
		// closest horizontal edge is the top one
		hseg_p1.y = hseg_p2.y = tr_edges.y;
		hseg_p1.x = ll_edges.x;
		hseg_p2.x = tr_edges.x;
		}
	// construct the two closest segments
	Line_segment vseg(vseg_p1, vseg_p2);
	Line_segment hseg(hseg_p1, hseg_p2);
	
	// compute the distances, and return whichever is least
	double dist_vseg = vseg.distance_from_segment(p);
	double dist_hseg = hseg.distance_from_segment(p);
	return (dist_vseg < dist_hseg) ? dist_vseg : dist_hseg;
}


/* This function uses a special-cased version of the "odd-even" rule for determining
whether a point is inside a shape. We intersect the shape with a ray 
from the point that is horizontal to the right.

We iterate through the line segments making up the side and count the number of 
intersections of a line segment and a ray from the point that is horizontal to the right. 
If the number is odd, the point is inside the shape. If zero or even, the 
point is outside the shape.  Special cases appear if the ray intersects one of the
vertices of the polygon.

In the real plane with computed coordinates, it is unlikely that
we will get a ray intersecting exactly on a vertex. But it could happen.
And of course, if the vertices and points have integral values, it is much more likely.
*/

double Polygon::distance_inside(const Point p) const
{
	const size_t num_pts = vertices.size();
	const size_t num_segments = num_pts - 1;
	double mindist = 0;
	int count = 0;
	double previous_dy = 0.0;

	// Extend a ray to the right of p and count the intersections
	// We go through the vertices a pair at a time, constructing each line segment.
	// i will be the "start" point p1, i + 1 the "stop" point, p2, of the segment.
	// The last point wraps around to use the first point.
	
	for (size_t ptno = 0; ptno < num_pts; ptno++) {
		// Construct the line segment - wrap around for the last one
		size_t second_ptno = (ptno == num_segments) ? 0 : ptno + 1;
		Line_segment ls(vertices[ptno], vertices[second_ptno]);
		// get the distance and keep track of the minimum
		// we need this information regardless of whether the point is inside or not,
		// and must always be the closest segment regardless.
		double dist = ls.distance_from_segment(p);
		if (ptno == 0 || dist < mindist) 
			mindist = dist;
		// skip if this segment is horizontal
		if(ls.is_horizontal())
			continue;
		// find intersection of horizontal ray with line segment
		double t = ls.parameter_given_y(p.y);
		// if ray doesn't intersect the line segment, skip it	
		if ( t < 0.0 || t > 1.0) {
			previous_dy = ls.get_dy();
			continue;
			}
		// get value of x for intersection point; skip if not to the right
		if (ls.x_given_parameter(t) < p.x) {
			previous_dy = ls.get_dy();
			continue;
			}

		// The ray intersects the line segment, but check for special cases.
		// If point of intersection is at p1, count only if 
		//    1. this is the first segment (previous dy not set)
		//    or 
		//    2. dy has the same sign as in the previous segment
		// is point of intersection at p1?
		if (t == 0.0) {
			if ((ptno == 0) || (ls.get_dy() * previous_dy > 0.0))
				count++;
			}
		// If point of intersection is at p2, do not count - will check at next p1;
		// otherwise, count the intersection and go to the next segment
		else if (t != 1.0)
			count++;
		previous_dy = ls.get_dy();
		}
	// p is inside if we counted an odd number of intersections.
	// return positive distance if p is inside, negative if outside
	if (count % 2)
		return mindist;
	else
		return -mindist;
}

// return the center of the bounding box
Point Polygon::get_center() const
{
	if(!cache_good)
		recompute();
	return center;
}

// return the Size of the bounding box
Size Polygon::get_size() const
{
	if(!cache_good)
		recompute();
	return size;
}

void Polygon::recompute() const
{
	if(vertices.empty()) {
		size.h = size.v = 0.;
		center.x = center.y = 0.;
		cache_good = true;
		return;
		}
	
	// find largest and smallest values in each dimension
	double xmin = vertices[0].x;
	double xmax = xmin;
	double ymin = vertices[0].y;
	double ymax = ymin;
	for(size_t i = 1; i < vertices.size(); ++i) {
		Point p = vertices[i];
		if(p.x < xmin)
			xmin = p.x;
		else if(p.x > xmax)
			xmax = p.x;
		if(p.y < ymin)
			ymin = p.y;
		else if(p.y > ymax)
			ymax = p.y;
		}
	size.h = fabs(xmax - xmin);
	size.v = fabs(ymax - ymin);
	center.x = xmin + size.h / 2.;
	center.y = ymin + size.v / 2.;
	return;
}


} // end namespace

------------------------------------------------------------
// Numeric_utilities.h
#ifndef NUMERIC_UTILITIES_H
#define NUMERIC_UTILITIES_H

#include <string>

// rounds x to an integer value returned in a double
double round_to_integer(double x);

// writes an integer value to a string
std::string int_to_string(int i);

// convert hours, minutes, and seconds to milliseconds (note long integer returned)
long time_convert(long hours, long minutes, long seconds);

// convert milliseconds into hours, minutes, seconds, and milliseconds
void time_convert(long time_ms, int& hours, int& minutes, int& seconds, int& milliseconds);

// convert milliseconds into hours, minutes, and double seconds
void time_convert(long time_ms, int& hours, int& minutes, double& seconds);

// compute the base 2 log of the supplied value
double logb2(double x);

double pitch_to_semitones(double pitch);



#endif

------------------------------------------------------------
// Numeric_utilities.cpp
#include "Numeric_utilities.h"

#include <string>
#include <sstream>
#include <cmath>

using std::string;
using std::ostringstream;
using std::log;

// return a double rounded to the nearest integer value
double round_to_integer(double x)
{
	int i = int(x + 0.5);
	return double(i);
}

string int_to_string(int i)
{
	ostringstream ss;
	ss << i;
	return ss.str();
}

// convert hours, minutes, and seconds to milliseconds (note long integer returned)
long time_convert(long hours, long minutes, long seconds)
{
	 return (seconds*1000 + minutes*60000 + hours*3600000);
}

// convert milliseconds into hours, minutes, seconds, and milliseconds
void time_convert(long time_ms, int& hours, int& minutes, int& seconds, int& milliseconds)
{
	hours = int(time_ms / 3600000);
	time_ms = time_ms % 3600000;
	minutes = int(time_ms / 60000);
	time_ms = time_ms % 60000;
	seconds = int(time_ms / 1000);
	time_ms = time_ms % 1000;
	milliseconds = int(time_ms);
}

// convert milliseconds into hours, minutes, and double seconds
void time_convert(long time_ms, int& hours, int& minutes, double& seconds)
{
	hours = int(time_ms / 3600000);
	time_ms = time_ms % 3600000;
	minutes = int(time_ms / 60000);
	time_ms = time_ms % 60000;
	int int_seconds = int(time_ms / 1000);
	time_ms = time_ms % 1000;
	int milliseconds = int(time_ms);
	seconds = double(int_seconds) + double(milliseconds) / 1000.;
}



// compute the base 2 log of the supplied value
double logb2(double x) 
{
	return (log(x) / log(2.0));
}


double pitch_to_semitones(double pitch)
{
    return logb2(pitch) * 12.;
}


------------------------------------------------------------
// PrdObsStatistics.h

//
//  Created by David Kieras on 4/17/18.
//

#ifndef PRDOBSSTATISTICS_H
#define PRDOBSSTATISTICS_H

#include "Program_constants.h"
#include "Statistics.h"
//#include "EPICLib/Statistics.h"

#include <iosfwd>

class PrdObsStatistics {
public:
    PrdObsStatistics();
    // read a file containing 3 X 8 RT values for CSF, COC, SHP, Each condition with values for Neg 3,6,12,18, Pos 3,5,12,18
    // followed by a 3 X 8 ER error rate values in the same order
    void load_obs_data(std::ifstream& infile, bool echo_print = true);

    void reset() {
        po_prop_errors.reset();
        po_rts.reset();
        }
    
    void update(Condition_e cond, Polarity_e polarity, int set_size, double pred_rt, double pred_err);

    void output(std::ostream& os) const;

private:
    double obsrts[3][2][4]; // condition X neg/pos X set size
    double obsproperrs[3][2][4]; // condition X neg/pos X set size
    
    PredObs_accumulator po_rts;
    PredObs_accumulator po_prop_errors;

    


};





#endif /* PRDOBSSTATISTICS_H */

------------------------------------------------------------
// PrdObsStatistics.cpp

//  SimpleCrowdingModel
//
//  Created by David Kieras on 4/17/18.
//

#include "PrdObsStatistics.h"
#include "Program_constants.h"
#include <string>
#include <sstream>
#include <cassert>
#include <iostream>
#include <fstream>


using namespace std;

PrdObsStatistics::PrdObsStatistics()
{
}

void PrdObsStatistics::load_obs_data(ifstream& infile, bool echo_print)
{
    for(int cond = 0; cond < 3; cond++) {
        for(int pol = 0; pol < 2; pol++) {
            for(int ssize = 0; ssize < 4; ssize++) {
                double rt;
                infile >> rt;
                assert(infile.good());
                obsrts[cond][pol][ssize] = rt;
                }
            }
        }
    for(int cond = 0; cond < 3; cond++) {
        for(int pol = 0; pol < 2; pol++) {
            for(int ssize = 0; ssize < 4; ssize++) {
                double er;
                infile >> er;
                assert(infile.good());
                obsproperrs[cond][pol][ssize] = er;
                }
            }
        }
    // echo print for check
    if(!echo_print)
    	return;
    for(int cond = 0; cond < 3; cond++) {
        for(int pol = 0; pol < 2; pol++) {
            for(int ssize = 0; ssize < 4; ssize++) {
                cout << obsrts[cond][pol][ssize] << ' ';
                cout << obsproperrs[cond][pol][ssize] << endl;
                }
            }
        }
}

void PrdObsStatistics::update(Condition_e cond, Polarity_e polarity, int set_size, double pred_rt, double pred_err)
{
    // set_size is 3/6/12/18 has to be transposed to 0/1/2/3
    int iset_size = 0; // value for 3
    if(set_size > 3)
        iset_size = set_size/6; // value for 6,12,18
    assert(iset_size >= 0 && iset_size <= 3);

    double obs_rt = obsrts[int(cond)][int(polarity)][iset_size];
    double obs_err = obsproperrs[int(cond)][int(polarity)][iset_size];
    
//    cout << int(cond) << ' ' << int(polarity) << ' ' << set_size << ' ' << iset_size << ' ' << obs_rt << ' ' << pred_rt << ' ' << obs_err << ' ' << pred_err << endl;
    // pred then observed because aare computed accordingly
    po_rts.update(pred_rt, obs_rt);
    po_prop_errors.update(pred_err, obs_err);
}

void PrdObsStatistics::output(ostream& os) const
{
    assert(po_rts.get_n() == po_prop_errors.get_n());
//    cout << po_rts.get_rsq() << "\t" << po_rts.get_avg_abs_rel_error() << "\t";
    os << po_rts.get_rsq() << "\t" << po_rts.get_avg_abs_error() << "\t" << po_rts.get_avg_abs_rel_error() << "\t";
    double rt_fom = po_rts.get_avg_abs_rel_error() / po_rts.get_rsq();
    os << rt_fom << "\t";
    os << po_prop_errors.get_rsq() << "\t" << po_prop_errors.get_avg_abs_error() << "\t" << po_prop_errors.get_avg_abs_rel_error() << "\t";
    double er_foma = po_prop_errors.get_avg_abs_error() / po_prop_errors.get_rsq();
    os <<  er_foma << "\t";
    double er_fomr = po_prop_errors.get_avg_abs_rel_error() / po_prop_errors.get_rsq();
    os <<  er_fomr << "\t";
    
    double wa_fomr = (RT_FoM_weight_c * rt_fom + (1. - RT_FoM_weight_c) * er_fomr);
    os << wa_fomr << "\t" << po_prop_errors.get_n();
}




------------------------------------------------------------
// Program_constants.h
//  SimpleCrowdingModel
//
//  Created by David Kieras on 3/12/18.
//
#ifndef PROGRAM_CONSTANTS_H
#define PROGRAM_CONSTANTS_H

//#include "EPICLib/Geometry.h"
#include "Geometry.h"
#include <array>
namespace GU = Geometry_Utilities;

// Global Types for experiment simulation
enum class Condition_e {CSF, COC, SHP};
enum class Polarity_e {NEGATIVE, POSITIVE};
enum class Run_mode_e {EXPERIMENT, GRID_SEARCH};

// availability constants
constexpr double common_intercept_c = 0;
constexpr double common_sd_c = 0.5;
// other visual constants
constexpr double bouma_fraction_c = 0.5;
// constexpr double bouma_fraction_c = 1.0; // for testing - massive crowding results
constexpr double fovea_radius_c = 1.0; // same as standard_fovea_radius in EPIC architecture


// these are declared as global variables rather than constants
// so that they can be changed in grid searches.
// Their initial values are defined as the constexpr with _c names.

extern double color_slope_g;
//extern double color_crowding_slope_g;
extern double color_crowding_probability_g;
extern double orientation_slope_g;
extern double orientation_crowding_probability_g;
extern double shape_slope_g;
extern double shape_crowding_slope_g;
extern double shape_crowding_probability_g;
extern double visual_delay_time_g;
extern int max_n_fixations_g;

extern int trial_g; // the current trial number, global for debugging output

//constexpr double color_slope_c = 0.0;
//constexpr double color_slope_c = 0.05;
//constexpr double color_slope_c = 0.055;
//constexpr double color_slope_c = 0.075;
//constexpr double color_slope_c = 0.08;
//constexpr double color_slope_c = 0.09;
//constexpr double color_slope_c = 0.095;
//constexpr double color_slope_c = 0.10;
constexpr double color_slope_c = 0.11; // * last setting as of 4/26/23
//constexpr double color_slope_c = 0.115;
//constexpr double color_slope_c = 0.125;
//constexpr double color_slope_c = 0.13;
//constexpr double color_slope_c = 0.135;
//constexpr double color_slope_c = 0.145;
//constexpr double color_slope_c = 0.15;
//constexpr double color_slope_c = 0.16;
//constexpr double color_slope_c = 0.17;
//constexpr double color_slope_c = 0.175;
//constexpr double color_slope_c = 0.18;
//constexpr double color_slope_c = 0.20;
//constexpr double color_slope_c = 0.30;
//constexpr double color_slope_c = 0.40; // for check against SHP
//constexpr double color_slope_c = 0.425; // for check against SHP

//constexpr double color_crowding_probability_c = 0.0;
//constexpr double color_crowding_probability_c = 0.01;
//constexpr double color_crowding_probability_c = 0.015; // * last setting as of 4/26/23
constexpr double color_crowding_probability_c = 0.025;
//constexpr double color_crowding_probability_c = 0.05;
//constexpr double color_crowding_probability_c = 0.075;
//constexpr double color_crowding_probability_c = 0.08;
//constexpr double color_crowding_probability_c = 0.10;
//constexpr double color_crowding_probability_c = 0.125;
//constexpr double color_crowding_probability_c = 0.15;
//constexpr double color_crowding_probability_c = 0.175;
//constexpr double color_crowding_probability_c = 0.2;
//constexpr double color_crowding_probability_c = 0.25;
//constexpr double color_crowding_probability_c = 0.30;
//constexpr double color_crowding_probability_c = 0.40;
//constexpr double color_crowding_probability_c = 0.50;
//constexpr double color_crowding_probability_c = 0.625;
//constexpr double color_crowding_probability_c = 0.75;
//constexpr double color_crowding_probability_c = 0.8;
//constexpr double color_crowding_probability_c = 1.0;

//constexpr double orientation_slope_c = 0.0;
//constexpr double orientation_slope_c = 0.025;
//constexpr double orientation_slope_c = 0.05;
//constexpr double orientation_slope_c = 0.075;
//constexpr double orientation_slope_c = 0.10;
//constexpr double orientation_slope_c = 0.11;
//constexpr double orientation_slope_c = 0.15;
//constexpr double orientation_slope_c = 0.17;
// constexpr double orientation_slope_c = 0.18; // * last setting as of 4/26/23
//constexpr double orientation_slope_c = 0.19;
constexpr double orientation_slope_c = 0.20;
//constexpr double orientation_slope_c = 0.225;
//constexpr double orientation_slope_c = 0.25;
//constexpr double orientation_slope_c = 0.275;
//constexpr double orientation_slope_c = 0.30;
//constexpr double orientation_slope_c = 0.315;
//constexpr double orientation_slope_c = 0.325;
//constexpr double orientation_slope_c = 0.35;
//constexpr double orientation_slope_c = 0.375;
//constexpr double orientation_slope_c = 0.4;
//constexpr double orientation_slope_c = 0.425;

//constexpr double orientation_crowding_probability_c = .0;
//constexpr double orientation_crowding_probability_c = .005;
//constexpr double orientation_crowding_probability_c = .01;
//constexpr double orientation_crowding_probability_c = .015;
constexpr double orientation_crowding_probability_c = .025; // * last setting as of 4/26/23
//constexpr double orientation_crowding_probability_c = .05;
//constexpr double orientation_crowding_probability_c = .075;
//constexpr double orientation_crowding_probability_c = .08;
//constexpr double orientation_crowding_probability_c = .10;
//constexpr double orientation_crowding_probability_c = .15;
//constexpr double orientation_crowding_probability_c = .20;
//constexpr double orientation_crowding_probability_c = .25;
//constexpr double orientation_crowding_probability_c = .30;
//constexpr double orientation_crowding_probability_c = .4;
//constexpr double orientation_crowding_probability_c = .45;
//constexpr double orientation_crowding_probability_c = .5;
//constexpr double orientation_crowding_probability_c = 1.0;

// constexpr double shape_slope_c = 0.0; // * last setting as of 4/26/23
//constexpr double shape_slope_c = 0.01;
//constexpr double shape_slope_c = 0.05;
//constexpr double shape_slope_c = 0.1;
//constexpr double shape_slope_c = 0.15;
//constexpr double shape_slope_c = 0.20;
//constexpr double shape_slope_c = 0.30;
//constexpr double shape_slope_c = 0.315;
//constexpr double shape_slope_c = 0.325;
//constexpr double shape_slope_c = 0.35;
//constexpr double shape_slope_c = 0.375;
//constexpr double shape_slope_c = 0.40;
//constexpr double shape_slope_c = 0.4125;
//constexpr double shape_slope_c = 0.42;
constexpr double shape_slope_c = 0.425;
//constexpr double shape_slope_c = 0.435;
//constexpr double shape_slope_c = 0.45;
//constexpr double shape_slope_c = 0.50;
//constexpr double shape_slope_c = 0.525;
//constexpr double shape_slope_c = 0.535;
//constexpr double shape_slope_c = 0.545;
//constexpr double shape_slope_c = 0.55;
//constexpr double shape_slope_c = 0.60;
//constexpr double shape_slope_c = 0.70;
//constexpr double shape_slope_c = 0.80;


constexpr double shape_crowding_slope_c = 0.0; // used in VM5+
//constexpr double shape_crowding_slope_c = 1.; // used in VM5+

// constexpr double shape_crowding_probability_c = 0.0;  // * last setting as of 4/26/23
//constexpr double shape_crowding_probability_c = 0.005;
//constexpr double shape_crowding_probability_c = 0.01;
//constexpr double shape_crowding_probability_c = 0.02;
//constexpr double shape_crowding_probability_c = 0.025;
//constexpr double shape_crowding_probability_c = 0.05;
constexpr double shape_crowding_probability_c = 0.075;
//constexpr double shape_crowding_probability_c = 0.1;
//constexpr double shape_crowding_probability_c = 0.125;
//constexpr double shape_crowding_probability_c = 0.135;
//constexpr double shape_crowding_probability_c = 0.15;
//constexpr double shape_crowding_probability_c = 0.18;
//constexpr double shape_crowding_probability_c = 0.20;
//constexpr double shape_crowding_probability_c = 0.25;
//constexpr double shape_crowding_probability_c = 0.4;
//constexpr double shape_crowding_probability_c = 0.5;
//constexpr double shape_crowding_probability_c = 0.8;
//constexpr double shape_crowding_probability_c = 1.0;

//constexpr double visual_delay_time_c = 24.; // original value
//constexpr double visual_delay_time_c = 25.;
//constexpr double visual_delay_time_c = 50.;
constexpr double visual_delay_time_c = 100;
//constexpr double visual_delay_time_c = 200;

// error rate parameters
//constexpr double negative_baseline_error_rate_c = .00;
//constexpr double negative_baseline_error_rate_c = .001;
//constexpr double negative_baseline_error_rate_c = .0015;
//constexpr double negative_baseline_error_rate_c = .004;
//constexpr double negative_baseline_error_rate_c = .0045;    // COC 1279
//constexpr double negative_baseline_error_rate_c = .005;
//constexpr double negative_baseline_error_rate_c = .006;
//constexpr double negative_baseline_error_rate_c = .007;
//constexpr double negative_baseline_error_rate_c = .009;
//constexpr double negative_baseline_error_rate_c = .01;
constexpr double negative_baseline_error_rate_c = .014; // average of allsubs negative
//constexpr double negative_baseline_error_rate_c = .02;  // CSF_1279
//constexpr double negative_baseline_error_rate_c = .021;  // CSF_1279 corrected 4/27/23
//constexpr double negative_baseline_error_rate_c = .023;
//constexpr double negative_baseline_error_rate_c = .029;  // COC_58
//constexpr double negative_baseline_error_rate_c = .03; // COC 458 * last setting as of 4/26/23
// constexpr double negative_baseline_error_rate_c = .032; // COC 458 * last setting as of 4/27/23

constexpr double positive_baseline_error_rate_c = negative_baseline_error_rate_c;
//constexpr double positive_baseline_error_rate_c = .005;
//constexpr double positive_baseline_error_rate_c = .006;
//constexpr double positive_baseline_error_rate_c = .007;
//constexpr double positive_baseline_error_rate_c = .01;
//constexpr double positive_baseline_error_rate_c = .011; // Avg ER+ for COC_1279 ss3,6
//constexpr double positive_baseline_error_rate_c = .014; // average of allsubs negative
//constexpr double positive_baseline_error_rate_c = .02;
//constexpr double positive_baseline_error_rate_c = .03;
//constexpr double positive_baseline_error_rate_c = .04;
//constexpr double positive_baseline_error_rate_c = .05;



// Strategy controls

// 5a: max_n_fixations_c = 99; confirm_positive_c = true; confirm_negative_c = true;
// almost 6b: max_n_fixations_c = comnfix; confirm_positive_c = true; confirm_negative_c = true;
// 7: max_n_fixations_c = 0; confirm_positive_c, confirm_negative_c, ignored but set to false;


//constexpr int max_n_fixations_c = 99; // effectively unlimited fixations // * last setting as of 4/26/23
//constexpr int max_n_fixations_c = 12;
//constexpr int max_n_fixations_c = 11;
//constexpr int max_n_fixations_c = 10;
//constexpr int max_n_fixations_c = 9;
//constexpr int max_n_fixations_c = 8;
//constexpr int max_n_fixations_c = 7;
//constexpr int max_n_fixations_c = 6;
//constexpr int max_n_fixations_c = 5;
//constexpr int max_n_fixations_c = 4;
//constexpr int max_n_fixations_c = 3;
//constexpr int max_n_fixations_c = 2;
//constexpr int max_n_fixations_c = 1;
constexpr int max_n_fixations_c = 0; // eyes kept on fixation point

// Probability that actual max number of fixations will be max_n_fixations_c+1
constexpr double max_n_fixations_mixture_probability_c = 0.0; // * last setting as of 4/26/23
//constexpr double max_n_fixations_mixture_probability_c = 0.25;
//constexpr double max_n_fixations_mixture_probability_c = 0.5;
//constexpr double max_n_fixations_mixture_probability_c = 0.75;
//constexpr double max_n_fixations_mixture_probability_c = 0.90;
//constexpr double max_n_fixations_mixture_probability_c = 0.95;

// calculate the average max_n_fixations value based on the _g value and the mixture probability
// assuming the mixture probability is for incrementing the base value by +1
double calculate_average_max_n_fixations();

constexpr bool confirm_positive_c = false;   // make confirmation eye move before positive response
constexpr bool confirm_negative_c = false;   // make confirmation eye move before negative response

constexpr double close_enough_to_bypass_confirm_positive_c = fovea_radius_c;  // default Strategy 8, 9 used this // * last setting as of 4/26/23
//constexpr double close_enough_to_bypass_confirm_positive_c = 2.0;
//constexpr double close_enough_to_bypass_confirm_positive_c = 3.0;
//constexpr double close_enough_to_bypass_confirm_positive_c = 0;  // almost never skip

constexpr int enough_fixations_to_bypass_confirm_negative_c = 1;   // default strategy 8. 9 used this // * last setting as of 4/26/23
//constexpr int enough_fixations_to_bypass_confirm_negative_c = 2;
//constexpr int enough_fixations_to_bypass_confirm_negative_c = 3;
//constexpr int enough_fixations_to_bypass_confirm_negative_c = 99; // never skip

constexpr int too_many_illusory_targets_c = 99; // use value of 99 to effectively turn off this check // * last setting as of 4/26/23
//constexpr int too_many_illusory_targets_c = 2;
//constexpr int too_many_illusory_targets_c = 3;
// constexpr int too_many_illusory_targets_c = 4;
//constexpr int too_many_illusory_targets_c = 5;

// all of below * last setting as of 4/26/23
// visual model mode parameters
// The following specifies policy for allowing properties from crowding objects to overwrite
// properties of an object that are in perceptual store state
enum class Perceptual_store_scramble_policy_e {NO_OVERWRITES, NON_BLANK_OVERWRITES, ALL_OVERWRITES};
constexpr auto perceptual_store_scramble_policy = Perceptual_store_scramble_policy_e::NON_BLANK_OVERWRITES;
constexpr bool apply_crowding_effect = true;
constexpr bool unitary_shape_c = true;
constexpr bool saccade_noise_present_c = true;

// below are initialized to false, and enabled by Strategy if it needs them
extern bool compute_Shape_information_gain_g;
extern bool compute_Feature_information_gain_g; // both Color and Orientation

// Observed data file
constexpr const char* const infile_name_c = "ObsData_AllSubjects.txt";
//constexpr const char* const infile_name_c = "OBsData_All_3610_37.txt";
//constexpr const char* const infile_name_c = "ObsData_48_310_129.txt";
//constexpr const char* const infile_name_c = "ObsData_356_1279_37.txt";
// constexpr const char* const infile_name_c = "ObsData_1279_458_458.txt";
//constexpr const char* const infile_name_c = "ObsData_1279_58_458.txt";


// Experiment control

constexpr auto conditions = {Condition_e::CSF, Condition_e::COC, Condition_e::SHP};
//constexpr auto conditions = {Condition_e::CSF, Condition_e::COC};
//constexpr auto conditions = {Condition_e::CSF};
//constexpr auto conditions = {Condition_e::COC};
//constexpr auto conditions = {Condition_e::SHP};

constexpr auto polarities = {Polarity_e::NEGATIVE, Polarity_e::POSITIVE};
//constexpr auto polarities = {Polarity_e::NEGATIVE};
//constexpr auto polarities = {Polarity_e::POSITIVE};

constexpr auto set_sizes = {3,6,12,18};
//constexpr auto set_sizes = {6};
//constexpr auto set_sizes = {12};
//constexpr auto set_sizes = {18};


// Simulation and Statistics
constexpr Run_mode_e run_mode_c = {Run_mode_e::EXPERIMENT};
//constexpr Run_mode_e run_mode_c = {Run_mode_e::GRID_SEARCH};
//constexpr bool run_experiment_c = false; // true for run experiment, false for do grid search
// number of trials/condition/polarity/setsize to run
//constexpr int n_trials_c = 3;
//constexpr int n_trials_c = 5;
//constexpr int n_trials_c = 10;
//constexpr int n_trials_c = 50;
//constexpr int n_trials_c = 100;
//constexpr int n_trials_c = 1000;
//constexpr int n_trials_c = 100000;
//constexpr int n_trials_c = 500000;
constexpr int n_trials_c = 1000000;

// weight of RT in calculating weighted average FoM for RT and ER [0, 1]
constexpr double RT_FoM_weight_c = 4./5.; // give equal weight for RT FoM ~ 5 and ER FoM ~20


// Additional visual model mode parameters
// SHP representation
// change below together
//constexpr int n_shape_parts_c = 7; // 7-segments
constexpr int n_shape_parts_c = 2;
using Shape_parts_t = std::array<char, n_shape_parts_c>;
//const std::array<int, n_shape_parts_c> part_indices_c{0,1,2,3,4,5,6}; // 7-segments
const std::array<int, n_shape_parts_c> part_indices_c{0,1};
// instead of having a separate variable, shape_slope_g is used when non-unitary shape parts are being detected
//const GU::Size shape_part_size_c = GU::Size(2.1, 0.3); // seven-segment size
const GU::Size shape_part_size_c = GU::Size(1.5, 1.35); // for 2-part shape: half of total shape height, same width




// debugging output
constexpr bool processed_visual_trace_output = false;
constexpr bool detailed_visual_trace_output = false;
constexpr bool crowding_visual_trace_output = false;
constexpr bool target_status_trace_output = false;
constexpr bool strategy_trace_output = false;
constexpr bool compute_and_output_display_statistics = false;

// time parameters
//constexpr double visual_encoding_time_c = 50.; // ms - not currently used - see visual_delay_time
//constexpr double same_key_motor_response_time_c = 50. + 25.;
//constexpr double different_key_motor_response_time_c = 100. + 50. + 25.;
// average of these to permit blocked trials
constexpr double same_key_motor_response_time_c = 50. + 50. + 25.;
constexpr double different_key_motor_response_time_c = 50. + 50. + 25.;
constexpr double em_time_intercept_c = 21.;
constexpr double em_time_slope_c = 2.2; // ms/dva
constexpr int rule_execution_time_c = 50; // ms/cycle


#endif



------------------------------------------------------------
// Program_constants.cpp

//  SimpleCrowdingModel
//
//  Created by David Kieras on 7/6/18.
//

#include "Program_constants.h"

// initialize these globals to the corresponding constant value
double color_slope_g = color_slope_c;
double color_crowding_probability_g = color_crowding_probability_c;
double orientation_slope_g = orientation_slope_c;
double orientation_crowding_probability_g = orientation_crowding_probability_c;
double shape_slope_g = shape_slope_c;
double shape_crowding_slope_g = shape_crowding_slope_c;
double shape_crowding_probability_g = shape_crowding_probability_c;
double visual_delay_time_g = visual_delay_time_c;
int max_n_fixations_g = max_n_fixations_c;

// below are enabled by Strategy if it needs them
bool compute_Shape_information_gain_g = false;
bool compute_Feature_information_gain_g = false; // both Color and Orientation

int trial_g = 0;

// calculate the average max_n_fixations value based on the _g value and the mixture probability
// assuming the mixture probability is for incrementing the base value by +1
double calculate_average_max_n_fixations()
{
    return (max_n_fixations_mixture_probability_c > 0.0) ?
        (max_n_fixations_g + max_n_fixations_mixture_probability_c) :
        max_n_fixations_g;
}



------------------------------------------------------------
// Random_utilities.h
//
//  Random_utilities.h
//
//  Created by David Kieras on 4/19/13.
//
//

#ifndef RANDOM_UTILITIES
#define RANDOM_UTILITIES

#include <random>
#include <vector>

// specify type of random engine here
typedef std::mt19937 Random_engine_t;
//typedef std::mt19937_64 Random_engine_t;
//typedef std::minstd_rand0 Random_engine_t;
//typedef std::knuth_b Random_engine_t;
//typedef std::default_random_engine Random_engine_t;

// accessor for the engine - same engine used for all randomization calls
Random_engine_t& get_Random_engine();

void set_random_number_generator_seed(unsigned long seed);
void reset_random_number_generator_seed();

/* Random variable generation */
// Returns a random integer in the range 0 ... range - 1 inclusive
int random_int(int range);
// Returns true with probability p
bool biased_coin_flip(double p);


double unit_uniform_random_variable();
// return a random variable that is uniformly distributed
// on each side of the mean +/- the deviation
double uniform_random_variable(double mean, double deviation);
double unit_normal_random_variable();
// do not call this function if sd == 0
double normal_random_variable(double mean, double sd);

// do not call this function if sd == 0
// return a normal variate in the range of +/- limit * sd,
// e.g. if limit == 3, then value is within +/- 3 sd of mean
double trimmed_normal_random_variable(double mean, double sd, double limit);

// return a normal random variable constrained to be > 0 - resample if not
double positive_normal_random_variable(double mean, double sd);

// this form uses a rate parameter, the reciprocal of the scale or survival parameter.
// the mean of this distribution is 1/rate
double exponential_random_variable_rate(double rate);

// this form uses a mean (expected survival time) parameter, the reciprocal of the rate parameter.
// the mean of this distribution is the expected survival time = 1/rate
double exponential_random_variable_mean(double mean);

//double floored_exponential_random_variable(double theta, double floor);
//double gamma_random_variable(double theta, int n);


// This form of log-normal uses the conventional mu and sigma as parameters,
// which govern location and spread of x, the random variable.
// x is constrained to be positive, with both low and high values less
// likely than values in the middle.
// mu is the mean of log(x), and sigma is the std. dev. of log(x).
// For random time values, make mu and sigma in seconds
// and multiply the result by 1000 to get milliseconds
double lognormal_random_variable_mu(double mu, double sigma);

// This form of log-normal uses median and spread as parameters,
// which govern location and spread of x, the random variable.
// x is constrained to be positive, with both low and high values less
// likely than values in the middle.
// The location parameter is the median of x, and spread is the std. dev. of log(x).
// Forbes, C., Evans, M. Hastings, N., Peacock, B.(2011).
// Statistical Distributions (4th ed). New York:Wiley. P. 131-134.
// For random time values, make mu and sigma in seconds
// and multiply the result by 1000 to get milliseconds
double lognormal_random_variable_median(double median, double spread);


// This form of Weibull distribution uses 3 prameters:
// location, scale (spread), and shape.
// Forbes, C., Evans, M. Hastings, N., Peacock, B.(2011).
// Statistical Distributions (4th ed). New York:Wiley. P. 193-201.
double weibull_random_variable(double location, double scale, double shape);


// This class generates a random value for an arbitrary discrete distribution described
// in the input vector of size n; it returns a value (0, n] using the probabilities in the vector.
// The values can be any probabilities, including zero, as long as they sum to one.
// The complex setup is done by the constructor and should not be repeated unnecessarily; getting a random value is then very fast.
class Discrete_distribution {
public:
    // copy the caller's vector because it will be modified during setup
    Discrete_distribution(std::vector<double> probabilities);
    
    // return a random index based on the probabilities
    int get_random_value() const;

private:
    int n;  // the number of probabilities supplied
    // these are the two tables used by the Alias Method
    std::vector<double> probability;
    std::vector<int> alias;

};

// Returns true with a probability = p
bool uniform_detection_function(double p);

// As x increases, the probability that true is returned increases according to a Normal dbn from 0. to 1.0.
bool gaussian_detection_function(double x, double mean, double sd);

// With lapse_probability, return false; else return the gaussian_detection_function result.
bool lapsed_gaussian_detection_function(double x, double mean, double sd, double lapse_probability);

// As x increases, the probability that true is returned increases according to a Normal dbn from base to 1.0.
bool based_gaussian_detection_function(double x, double base, double mean, double sd);

// As x increases, the probability that true is returned increases according to a Normal dbn from 0 to cap.
bool capped_gaussian_detection_function(double x, double cap, double mean, double sd);

// As x increases, the probability that true is returned increases according to an exponential dbn from 0 to 1.0
// the lambda parameter is assumed to be the "small" definition - 
// e.g. a value of lambda = .5 gives p(x < 1) = .3868
// cf lambda = 2 giving p(x < 1) = .8706
bool exponential_detection_function(double x, double lambda);

// As x increases, the probability that true is returned increases according to a exponential dbn from base to 1.0.
bool based_exponential_detection_function(double x, double base, double lambda);

// given two z-scores, return the cumulative probability from the bivariate normal.
// e.g. (-3, +3) returns  0.001348; (+3, +3) returns 0.997302.
// (+1, +2) and (+2, +1) return 0.822204

double get_bivariate_normal_cdf(double z1, double z2);


#endif


------------------------------------------------------------
// Random_utilities.cpp
//
//  Random_utilities.cpp
//
//  Created by David Kieras on 4/19/13.
//
//

#include "Random_utilities.h"
//#include "Assert_throw.h"
#include <random>
#define USE_ASSERT_THROW

#ifdef USE_ASSERT_THROW
#include "Assert_throw.h"
#define ASSERT Assert
#endif

#ifdef USE_CASSERT
#include <cassert>
#define ASSERT assert
#endif


using namespace std;

// This code uses C++11's <random> library
// functions declared here for inline purposes

/*
typedef std::mt19937 Random_engine_t;
//typedef std::minstd_rand0 Random_engine_t;
//typedef std::knuth_b Random_engine_t;
//typedef std::default_random_engine Random_engine_t;

// The static local variable in this function is the global random engine, hidden inside the unnamed namespace
namespace {
Random_engine_t& get_Random_engine()
{
    static Random_engine_t engine;
    return engine;
}
}
*/


Random_engine_t& get_Random_engine()
{
    static Random_engine_t engine;
    return engine;
}


void set_random_number_generator_seed(unsigned long seed_)
{
    get_Random_engine().seed(static_cast<Random_engine_t::result_type>(seed_));
}

void reset_random_number_generator_seed()
{
    get_Random_engine().seed();
}

/* Random variable generation - using C++11 <random> */


// this returns a random integer in the range [0, range)
// that is, 0 ... range - 1

int random_int(int range)
{
    // uniform_int_distribution returns a <= x <= b
    static std::uniform_int_distribution<int> uid;
    using parm_t = decltype(uid)::param_type;
    int result = uid(get_Random_engine(), parm_t(0, range - 1));
    ASSERT(result != range);
    return result;
}

// Returns true with probability p
bool biased_coin_flip(double p)
{
	double rv = unit_uniform_random_variable();
	return (rv <= p);
}

double unit_uniform_random_variable()
{
    static std::uniform_real_distribution<double> uurd(0., 1.);
    return uurd(get_Random_engine());
}

// return a random variable that is uniformly distributed
// on each side of the mean +/- the deviation
double uniform_random_variable(double mean, double deviation)
{
	return 2. * deviation * unit_uniform_random_variable() - deviation + mean;
}

double unit_normal_random_variable()
{
    static std::normal_distribution<double> und(0., 1.);
    return und(get_Random_engine());
}


double normal_random_variable(double mean, double sd)
{
    static std::normal_distribution<double> nd;
    using parm_t = decltype(nd)::param_type;
    return nd(get_Random_engine(), parm_t(mean, sd));
}

// do not call this function if sd == 0
// return a normal variate in the range of +/- limit * sd,
// e.g. if limit == 3, then value is within +/- 3 sd of mean
double trimmed_normal_random_variable(double mean, double sd, double limit)
{
    static std::normal_distribution<double> nd;
    using parm_t = decltype(nd)::param_type;
	double x, deviation;
	do {
		x = nd(get_Random_engine(), parm_t(mean, sd));
		deviation = abs(x - mean)/sd;
		} while (deviation > limit);
	return x;
}

// return a normal random variable constrained to be > 0 - resample if not
double positive_normal_random_variable(double mean, double sd)
{
    static std::normal_distribution<double> nd;
    using parm_t = decltype(nd)::param_type;
	double x;
	do {
		x = nd(get_Random_engine(), parm_t(mean, sd));
		} while (x <= 0.0);
	return x;
}


double exponential_random_variable_rate(double rate)
{
// std::exponential_distribution uses a rate parameter = 1/mean
 //	return -theta * log(unit_uniform_random_variable());
    static std::exponential_distribution<double> ed;
    using parm_t = decltype(ed)::param_type;
    return ed(get_Random_engine(), parm_t(rate));
}

double exponential_random_variable_mean(double mean)
{
// std::exponential_distribution uses a rate parameter = 1/mean
 //return -mean * log(unit_uniform_random_variable());
    static std::exponential_distribution<double> ed;
    using parm_t = decltype(ed)::param_type;
    return ed(get_Random_engine(), parm_t(1./mean));
}



// This form of log-normal uses the conventional mu and sigma as parameters,
// which govern location and spread of x, the random variable.
// x is constrained to be positive, with both low and high values less
// likely than values in the middle.
// mu is the mean of log(x), and sigma is the std. dev. of log(x).
// For random time values, make mu and sigma in seconds
// and multiply the result by 1000 to get milliseconds
double lognormal_random_variable_mu(double mu, double sigma)
{
    // below uses stdlib parameterization
    static std::lognormal_distribution<double> lnd;
    using parm_t = decltype(lnd)::param_type;
    return lnd(get_Random_engine(), parm_t(mu, sigma));
    
}

// This form of log-normal uses median and spread as parameters,
// which govern location and spread of x, the random variable.
// x is constrained to be positive, with both low and high values less
// likely than values in the middle.
// The location parameter is the median of x, and spread is the std. dev. of log(x).
// Forbes, C., Evans, M. Hastings, N., Peacock, B.(2011).
// Statistical Distributions (4th ed). New York:Wiley. P. 131-134.
// For random time values, make mu and sigma in seconds
// and multiply the result by 1000 to get milliseconds
double lognormal_random_variable_median(double median, double spread)
{
// See Forbes, C., Evans, M. Hastings, N., Peacock, B.(2011).
// Statistical Distributions (4th ed). New York:Wiley. P. 131-134.
	return median * exp(spread * unit_normal_random_variable());
}

// This form of Weibull distribution uses 3 prameters:
// location, scale (spread), and shape.
// Forbes, C., Evans, M. Hastings, N., Peacock, B.(2011).
// Statistical Distributions (4th ed). New York:Wiley. P. 193-201.
double weibull_random_variable(double location, double scale, double shape)
{
	double base = -log(unit_uniform_random_variable());
	double exponent = 1. / shape;
	double variate = pow(base, exponent);
	double x = location + scale * variate;
	return x;
}

// This code is based on the alias method as described by Keith Schwarz at
// http://keithschwarz.com/darts-dice-coins/
// and is very closely directly transcribed into C++ from his Java code at
// http://www.keithschwarz.com/interesting/code/?dir=alias-method

// caller's vector needs to be copied
Discrete_distribution::Discrete_distribution(vector<double> probabilities) :
    // initialize the two arrays to have the same number of cells as the input
    n(int(probabilities.size())), probability(n), alias(n)
{
    ASSERT(n);   // must be non-empty
    
    // average probability
    double average = 1.0 / n;
    
    // sort the probabilites into a large and small list
    // they might be of different lengths
    vector<int> small;
    vector<int> large;
    
    for(int i = 0; i < n; i++) {
        if(probabilities[i] >= average)
            large.push_back(i);
        else
            small.push_back(i);
        }
    
    // go through the lists together as long as there are entries in both,
    // and get the index of a large and small probability
    while(!small.empty() && !large.empty()) {
        int less = small.back();
        small.pop_back();
        int more = large.back();
        large.pop_back();
        
        // scale them up so that 1/n has weight 1.0
        probability[less] = probabilities[less] * n;
        alias[less] = more;
        // decrease the larger probability by the appropriate amount
        probabilities[more] = probabilities[more] + probabilities[less] - average;
        
        // put the new probability in the large or small list and continue the loop.
        if(probabilities[more] >= 1.0 / n)
            large.push_back(more);
        else
            small.push_back(more);
        }
    
    // at this point there may be left-over values in either list, so set those
    // indices to 1 in the probability table.
    while(!small.empty()) {
        probability[small.back()] = 1.0;
        small.pop_back();
        }
    while(!large.empty()) {
        probability[large.back()] = 1.0;
        large.pop_back();
        }
}

// return a random integer (0, n) based on the original probabilities
// by first picking a "column" in the 2D table, and then using probability
// in that "column" to flip a biased coin, and return either the "column" index
// or the alias at that index.
int Discrete_distribution::get_random_value() const
{
    int column = random_int(n);
    if(biased_coin_flip(probability[column]))
        return column;
    else
        return alias[column];
}


// returns true with a probability = p
bool uniform_detection_function(double p)
{
	double rv = unit_uniform_random_variable();
	return (rv <= p);
}

// As x increases, the probability that true is returned increases according to a Normal dbn from 0. to 1.0.
bool gaussian_detection_function(double x, double mean, double sd)
{
	double z = (x - mean) / sd;	// compute z-score
	double rv = unit_normal_random_variable();
	return (rv <= z);
}

// With lapse_probability, return false; else return the gaussian_detection_function result.
bool lapsed_gaussian_detection_function(double x, double mean, double sd, double lapse_probability)
{
    if(biased_coin_flip(lapse_probability))
        return false;
    return gaussian_detection_function(x, mean, sd);
}


// As x increases, the probability that true is returned increases according to a Normal dbn from base to 1.0.
bool based_gaussian_detection_function(double x, double base, double mean, double sd)
{
	double rv = unit_uniform_random_variable();
	if(rv <= base)
		return true;
	return gaussian_detection_function(x, mean, sd);
}

// As x increases, the probability that true is returned increases according to a Normal dbn from 0 to cap.
bool capped_gaussian_detection_function(double x, double cap, double mean, double sd)
{
	double rv = unit_uniform_random_variable();
	if(rv <= 1.0 - cap)
		return false;
	return gaussian_detection_function(x, mean, sd);
}

// As x increases, the probability that true is returned increases according to an exponential dbn from 0 to 1.0
// the lambda parameter is assumed to be the "small" definition - 
// e.g. a value of lambda = .5 gives p(x < 1) = .3868
// cf lambda = 2 giving p(x < 1) = .8706
bool exponential_detection_function(double x, double lambda)
{
	double rv = -log(unit_uniform_random_variable()) / lambda;
	return (rv <= x);
}

// As x increases, the probability that true is returned increases according to a exponential dbn from base to 1.0.
bool based_exponential_detection_function(double x, double base, double lambda)
{
	double rv = unit_uniform_random_variable();
	if(rv <= base)
		return true;
	return exponential_detection_function(x, lambda);
}



// given two z-scores, return the cumulative probability from the bivariate normal.
// e.g. (-3, +3) returns  0.001348; (+3, +3) returns 0.997302.
// (+1, +2) and (+2, +1) return 0.822204

/* implementation for get_bivariate_normal_cdf as fast table look-up */
// these constants should all have internal linkage

const int table_size = 25;
const double table[table_size][table_size] = {
/*          0         1         2         3         4         5         6         7         8         9         10        11        12        13        14        15        16        17        18        19        20        21        22        23        24 */
/* 0  */{0.000002, 0.000004, 0.000008, 0.000017, 0.000031, 0.000054, 0.000090, 0.000143, 0.000214, 0.000306, 0.000416, 0.000542, 0.000675, 0.000808, 0.000933, 0.001044, 0.001136, 0.001207, 0.001260, 0.001296, 0.001319, 0.001333, 0.001342, 0.001346, 0.001348}, 
/* 1  */{0.000004, 0.000009, 0.000019, 0.000036, 0.000068, 0.000119, 0.000199, 0.000315, 0.000473, 0.000675, 0.000919, 0.001196, 0.001490, 0.001784, 0.002060, 0.002304, 0.002507, 0.002665, 0.002781, 0.002860, 0.002912, 0.002943, 0.002961, 0.002971, 0.002976}, 
/* 2  */{0.000008, 0.000019, 0.000039, 0.000076, 0.000141, 0.000249, 0.000415, 0.000656, 0.000985, 0.001407, 0.001916, 0.002492, 0.003105, 0.003718, 0.004294, 0.004802, 0.005224, 0.005554, 0.005795, 0.005961, 0.006068, 0.006134, 0.006171, 0.006191, 0.006201}, 
/* 3  */{0.000017, 0.000036, 0.000076, 0.000149, 0.000278, 0.000490, 0.000817, 0.001292, 0.001939, 0.002770, 0.003772, 0.004906, 0.006112, 0.007319, 0.008453, 0.009454, 0.010285, 0.010933, 0.011408, 0.011735, 0.011946, 0.012075, 0.012149, 0.012188, 0.012208}, 
/* 4  */{0.000031, 0.000068, 0.000141, 0.000278, 0.000518, 0.000911, 0.001520, 0.002404, 0.003609, 0.005156, 0.007019, 0.009129, 0.011375, 0.013621, 0.015731, 0.017594, 0.019141, 0.020347, 0.021230, 0.021839, 0.022233, 0.022472, 0.022609, 0.022682, 0.022719}, 
/* 5  */{0.000054, 0.000119, 0.000249, 0.000490, 0.000911, 0.001605, 0.002676, 0.004232, 0.006356, 0.009079, 0.012360, 0.016075, 0.020030, 0.023984, 0.027699, 0.030981, 0.033704, 0.035827, 0.037383, 0.038454, 0.039148, 0.039569, 0.039810, 0.039940, 0.040005}, 
/* 6  */{0.000090, 0.000199, 0.000415, 0.000817, 0.001520, 0.002676, 0.004463, 0.007058, 0.010599, 0.015140, 0.020613, 0.026809, 0.033404, 0.039998, 0.046195, 0.051667, 0.056208, 0.059749, 0.062344, 0.064131, 0.065287, 0.065991, 0.066392, 0.066608, 0.066717}, 
/* 7  */{0.000143, 0.000315, 0.000656, 0.001292, 0.002404, 0.004232, 0.007058, 0.011162, 0.016762, 0.023943, 0.032597, 0.042397, 0.052825, 0.063253, 0.073053, 0.081707, 0.088888, 0.094488, 0.098592, 0.101418, 0.103246, 0.104358, 0.104994, 0.105335, 0.105507}, 
/* 8  */{0.000214, 0.000473, 0.000985, 0.001939, 0.003609, 0.006356, 0.010599, 0.016762, 0.025171, 0.035956, 0.048951, 0.063667, 0.079328, 0.094988, 0.109704, 0.122700, 0.133484, 0.141893, 0.148056, 0.152300, 0.155046, 0.156716, 0.157670, 0.158182, 0.158441}, 
/* 9  */{0.000306, 0.000675, 0.001407, 0.002770, 0.005156, 0.009079, 0.015140, 0.023943, 0.035956, 0.051360, 0.069923, 0.090944, 0.113314, 0.135683, 0.156704, 0.175267, 0.190672, 0.202684, 0.211487, 0.217549, 0.221472, 0.223857, 0.225220, 0.225952, 0.226321}, 
/* 10 */{0.000416, 0.000919, 0.001916, 0.003772, 0.007019, 0.012360, 0.020613, 0.032597, 0.048951, 0.069923, 0.095195, 0.123814, 0.154269, 0.184723, 0.213342, 0.238614, 0.259586, 0.275941, 0.287925, 0.296178, 0.301518, 0.304766, 0.306622, 0.307618, 0.308121}, 
/* 11 */{0.000542, 0.001196, 0.002492, 0.004906, 0.009129, 0.016075, 0.026809, 0.042397, 0.063667, 0.090944, 0.123814, 0.161037, 0.200647, 0.240257, 0.277480, 0.310350, 0.337626, 0.358897, 0.374484, 0.385218, 0.392164, 0.396388, 0.398802, 0.400098, 0.400752}, 
/* 12 */{0.000675, 0.001490, 0.003105, 0.006112, 0.011375, 0.020030, 0.033404, 0.052825, 0.079328, 0.113314, 0.154269, 0.200647, 0.250000, 0.299353, 0.345731, 0.386686, 0.420672, 0.447175, 0.466596, 0.479970, 0.488625, 0.493888, 0.496895, 0.498510, 0.499325}, 
/* 13 */{0.000808, 0.001784, 0.003718, 0.007319, 0.013621, 0.023984, 0.039998, 0.063253, 0.094988, 0.135683, 0.184723, 0.240257, 0.299353, 0.358449, 0.413983, 0.463023, 0.503718, 0.535453, 0.558708, 0.574723, 0.585086, 0.591387, 0.594989, 0.596922, 0.597898}, 
/* 14 */{0.000933, 0.002060, 0.004294, 0.008453, 0.015731, 0.027699, 0.046195, 0.073053, 0.109704, 0.156704, 0.213342, 0.277480, 0.345731, 0.413983, 0.478120, 0.534758, 0.581758, 0.618410, 0.645268, 0.663763, 0.675732, 0.683010, 0.687169, 0.689402, 0.690529}, 
/* 15 */{0.001044, 0.002304, 0.004802, 0.009454, 0.017594, 0.030981, 0.051667, 0.081707, 0.122700, 0.175267, 0.238614, 0.310350, 0.386686, 0.463023, 0.534758, 0.598105, 0.650673, 0.691666, 0.721706, 0.742392, 0.755778, 0.763919, 0.768570, 0.771068, 0.772329}, 
/* 16 */{0.001136, 0.002507, 0.005224, 0.010285, 0.019141, 0.033704, 0.056208, 0.088888, 0.133484, 0.190672, 0.259586, 0.337626, 0.420672, 0.503718, 0.581758, 0.650673, 0.707861, 0.752457, 0.785137, 0.807641, 0.822204, 0.831060, 0.836120, 0.838838, 0.840209}, 
/* 17 */{0.001207, 0.002665, 0.005554, 0.010933, 0.020347, 0.035827, 0.059749, 0.094488, 0.141893, 0.202684, 0.275941, 0.358897, 0.447175, 0.535453, 0.618410, 0.691666, 0.752457, 0.799862, 0.834601, 0.858523, 0.874004, 0.883417, 0.888797, 0.891685, 0.893143}, 
/* 18 */{0.001260, 0.002781, 0.005795, 0.011408, 0.021230, 0.037383, 0.062344, 0.098592, 0.148056, 0.211487, 0.287925, 0.374484, 0.466596, 0.558708, 0.645268, 0.721706, 0.785137, 0.834601, 0.870849, 0.895810, 0.911963, 0.921785, 0.927398, 0.930412, 0.931933}, 
/* 19 */{0.001296, 0.002860, 0.005961, 0.011735, 0.021839, 0.038454, 0.064131, 0.101418, 0.152300, 0.217549, 0.296178, 0.385218, 0.479970, 0.574723, 0.663763, 0.742392, 0.807641, 0.858523, 0.895810, 0.921486, 0.938102, 0.948206, 0.953980, 0.957080, 0.958645}, 
/* 20 */{0.001319, 0.002912, 0.006068, 0.011946, 0.022233, 0.039148, 0.065287, 0.103246, 0.155046, 0.221472, 0.301518, 0.392164, 0.488625, 0.585086, 0.675732, 0.755778, 0.822204, 0.874004, 0.911963, 0.938102, 0.955017, 0.965304, 0.971181, 0.974338, 0.975931}, 
/* 21 */{0.001333, 0.002943, 0.006134, 0.012075, 0.022472, 0.039569, 0.065991, 0.104358, 0.156716, 0.223857, 0.304766, 0.396388, 0.493888, 0.591387, 0.683010, 0.763919, 0.831060, 0.883417, 0.921785, 0.948206, 0.965304, 0.975700, 0.981642, 0.984832, 0.986442}, 
/* 22 */{0.001342, 0.002961, 0.006171, 0.012149, 0.022609, 0.039810, 0.066392, 0.104994, 0.157670, 0.225220, 0.306622, 0.398802, 0.496895, 0.594989, 0.687169, 0.768570, 0.836120, 0.888797, 0.927398, 0.953980, 0.971181, 0.981642, 0.987619, 0.990829, 0.992449}, 
/* 23 */{0.001346, 0.002971, 0.006191, 0.012188, 0.022682, 0.039940, 0.066608, 0.105335, 0.158182, 0.225952, 0.307618, 0.400098, 0.498510, 0.596922, 0.689402, 0.771068, 0.838838, 0.891685, 0.930412, 0.957080, 0.974338, 0.984832, 0.990829, 0.994049, 0.995674}, 
/* 24 */{0.001348, 0.002976, 0.006201, 0.012208, 0.022719, 0.040005, 0.066717, 0.105507, 0.158441, 0.226321, 0.308121, 0.400752, 0.499325, 0.597898, 0.690529, 0.772329, 0.840209, 0.893143, 0.931933, 0.958645, 0.975931, 0.986442, 0.992449, 0.995674, 0.997302} 
/*        -3        -2.75     -2.5      -2.25     -2.0      -1.75     -1.5      -1.25     -1.0      -.75      -.50      -.25      0         +.25      +.5       +.75      +1.0      +1.25     +1.5      + 1.75    +2.0      +2.25      +2.50    +2.75     +3.0 */
};

// range of table is -3 to +3 z-score in each dimension, 25 cells, .25 in each cell
const double min_z = -3.;
const double max_z = +3;
const double cell_size = 0.25;

static inline int compute_index(double z)
{
	int i;
	if(z <= min_z)
		i = 0;
	else if (z >= max_z)
		i = table_size - 1;
	else
		i = int((z - min_z) / cell_size);
	ASSERT(i >= 0 && i < table_size);
	return i;
}

double get_bivariate_normal_cdf(double z1, double z2)
{
	int i1 = compute_index(z1);
	int i2 = compute_index(z2);
		
//	cout << "i1, i2: " << i1 << ' ' << i2 << endl;

	return table[i1][i2];

}

/*
int random_int(int range)
{
	int result;
	do {
		double x = double(rand())/RAND_MAX;
		// note! it is possible for x to end up as 1.00000...
		result = int(double(range) * x);
		// so resample if it happens
		} while (result == range);
		
//	Assert(result != range);
	return result;
//	return int(double(range) * x);
}
*/
/*
double unit_uniform_random_variable()
{
	double x;
	do {
		x = double(rand())/RAND_MAX;
		} while (x == 0.0 || x == 1.0);
	return x;
}
*/
/*
double unit_normal_random_variable()
{
	double sum = -6.0;
	for(int i = 0; i < 12; i++)
		sum += unit_uniform_random_variable();
	return sum;
}
*/
/*
double normal_random_variable(double mean, double sd)
{
	return mean + sd * unit_normal_random_variable();
}
*/
/*
;return a random variable that follows the exponential distribution
(DEFUN EXPONENTIAL-RANDOM-VARIABLE (THETA)
(* (- THETA) (LOG (UNIT-UNIFORM-RANDOM-VARIABLE)))
)

;Since EPIC times are measured in milliseconds, this delivers an integer
;Theta should be in millisecond units also!
(DEFUN EXPONENTIAL-RANDOM-INTEGER (THETA)
(TRUNCATE (EXPONENTIAL-RANDOM-VARIABLE THETA))
)

;return a random variable that follows the exponential distribution
;plus a minimum floor
(DEFUN FLOORED-EXPONENTIAL-RANDOM-VARIABLE (THETA FLOOR)
(+ (* (- THETA) (LOG (UNIT-UNIFORM-RANDOM-VARIABLE))) FLOOR)
)
*/

/*
double floored_exponential_random_variable(double theta, double floor)
{
	return floor + (-theta) * log(unit_uniform_random_variable());
}
*/

/*
;Since EPIC times are measured in milliseconds, this delivers an integer
;Theta and floor should be in millisecond units also!
(DEFUN FLOORED-EXPONENTIAL-RANDOM-INTEGER (THETA FLOOR)
(TRUNCATE (FLOORED-EXPONENTIAL-RANDOM-VARIABLE THETA FLOOR))
)

;return a gamma random variable with parameters N and THETA
(DEFUN GAMMA-RANDOM-VARIABLE(THETA N)
(LET ((SUM 0.))
(DOTIMES (I N)
  (INCF SUM (EXPONENTIAL-RANDOM-VARIABLE THETA)))
SUM
))
*/

/*
double gamma_random_variable(double theta, int n)
{
	double sum = 0.;
	for(int i = 0; i < n; i++)
		sum += exponential_random_variable(theta);
	return sum;
}
*/

/*
;Since EPIC times are measured in milliseconds, this delivers an integer
;Theta should be in millisecond units also!
(DEFUN GAMMA-RANDOM-INTEGER (THETA N)
(TRUNCATE (GAMMA-RANDOM-VARIABLE THETA N))
)

;return a log-normal random variable with median M and parameter S
(DEFUN LOG-NORMAL-RANDOM-VARIABLE(M S)
(* M (EXP (* S (UNIT-NORMAL-RANDOM-VARIABLE))))
)
*/
/*
;Since EPIC times are measured in milliseconds, this delivers an integer
;M and S should be in millisecond units also!
(DEFUN LOG-NORMAL-RANDOM-INTEGER (M S)
(TRUNCATE (LOG-NORMAL-RANDOM-VARIABLE M S))
)


;**** RANDOM VARIABLE GENERATORS ***
;this function returns a number in the open interval (0, 1)
;by excluding the zero possibility from (random n)
(DEFUN UNIT-UNIFORM-RANDOM-VARIABLE ()
(LET (X)
(LOOP
  (SETQ X (RANDOM 100000))
  (IF (> X 0) (RETURN)))
(/ (FLOAT X) 100000.)
))

(DEFUN RANDOM-UNIFORM-DEVIATION (R)
(- (* R (UNIT-UNIFORM-RANDOM-VARIABLE)) (/ R 2.))
)

;return a unit normal variate
(DEFUN UNIT-NORMAL-RANDOM-VARIABLE ()
(LET ((SUM 0.))
(DOTIMES (I 12)
  (INCF SUM (UNIT-UNIFORM-RANDOM-VARIABLE)))
(- SUM 6.0)
))

;collect n samples by evaluating dbn-expression, and count them
;in the bins defined by max-value and number-bins.  All values
;above the top bin are included in the top bin.
;assumes that the distribution has a theoretical minimum value of zero, so
;only positive values will be obtained.
;sample:
;  (test-distribution '(LOG-NORMAL-RANDOM-INTEGER 3500 .25) 5000 20 1000) ;note quote!
(defun test-distribution (dbn-expression max-value number-bins n)
(let (bins bin-index step-size x min_x max_x (p 0.) (sum 0.) (totalp 0.))

(setq bins (make-array number-bins :initial-element 0))
(setq step-size (/ max-value number-bins))
(dotimes (i n)
  (setq x (eval dbn-expression))
  (incf sum x)
  (if (or (not min_x) (< x min_x)) (setq min_x x))
  (if (or (not max_x) (> x max_x)) (setq max_x x))
  (setq bin-index (truncate (/ x step-size)))
  (if (> bin-index (1- number-bins)) (setq bin-index (1- number-bins)))
  (incf (aref bins bin-index))
)

(format t "     x    f(x)   F(x)~%")
(setq x step-size)
(dotimes (i number-bins)
  (setq p (/ (aref bins i) n))
  (incf totalp p)
  (format t "   ~4,1f  ~5,3f  ~5,3f~%" x p totalp)
  (setq x (+ x step-size))
)
(format t "~%mean = ~5,2f, min = ~5,2f, max = ~5,2f~%" 
        (/ sum n) min_x max_x)
))

*/


------------------------------------------------------------
// Search_object.h
//  
//
//  Created by David Kieras on 3/7/18.
//

#ifndef SEARCH_OBJECT_H
#define SEARCH_OBJECT_H

//#include "EPICLib/Geometry.h"
#include "Geometry.h"
#include "Program_constants.h"
#include <iosfwd>
#include <vector>
#include <array>
#include <map>

namespace GU = Geometry_Utilities;

struct Search_object {
    const char* name{"\0"}; // object type string, distractor vs target
    int ID{-1};    // physical object ID
    GU::Point location{};
    GU::Size size{};
    double eccentricity{-1};
    double nearest_neighbor_distance = -1; // distance to nearest neighbor
    double critical_spacing = -1.; // current critical spacing for finding crowding objects
    std::vector<int> crowder_group; // indices into search_objects of crowding objects; always includes this object
//    enum class Crowding_state_e {Unknown, Crowded, Uncrowded};
//    Crowding_state_e crowding_state = Crowding_state_e::Unknown;
	// below not used in Visual_modelV2d,e
    bool current_crowding = false; // true if this object either crowds, or is crowded by another object.
    bool known_perceptual_properties = false; // true if this object's properties have become known while uncrowded
    // below added circa 4/26/19 to go with VM2b
    double Shape_information_gain{0.0}; // the score for gain of information about shape if this object is fixated
    double Feature_information_gain{0.0}; // the score for gain of information about color and shape if this object is fixated

	// these are the true physical properties
	char phys_color{'.'};
    char phys_orientation{'.'};
    char phys_shape{'.'};

    bool is_target = false; // true if this is actually a target object - for ease in monitoring scrambling results

    // these are the sensory properties - change with eye movements
    char sens_color{'.'};
    char sens_orientation{'.'};
    char sens_shape{'.'};
 
    // these are the perceptual stored properties - subject to crowding, sticky past eye movements, but updated
    char perc_color{'.'};
    char perc_orientation{'.'};
    char perc_shape{'.'};
 
    Shape_parts_t sens_spat_config{spat_configNull};

   /* "seven segment" configuration
   arramge with horizontal first, top to bottom, then vertical top left,bottom left, top right, bottom right
        ---0---     target 2 is 0,1,2,b,4,5,b in numerical segment index order, 0,5,1,4,2 in stroke order
        |     |     distractor 5 is 0,1,2,3,b,b,6 in numerical segment index order, 0,3,1,6,2 in stroke order
        3     5
        |     |     0 - 6 are locations (i.e. centers of segments - left/center/right, top/middle/bottom, horizontal/vertical line segments
        |--1--|     both 2 and 5 have H top, middle, bottom, 2 also has V left top, V right bottom, 5 also has V right top, V left bottom
        |     |
        4     6     in array: 2: H,H,H,b,V,V,b  5: H,H,H,V,b,b,V
        |     |
         --2--
         assume line segments are .3 dva in thickness, characters are 1.5 X 2.7, so each segment is about .3 X 2.1 -> 1.2 dva avg size
    */
/*
    Shape_parts_t perc_spat_config{spat_configNull};    // "seven segment" spatial configuration:
    static constexpr Shape_parts_t spat_configNull = {'.','.','.','.','.','.','.'};
    static constexpr Shape_parts_t spat_config2 = {'H','H','H','b','V','V','b'};
    static constexpr Shape_parts_t spat_config5 = {'H','H','H','V','b','b','V'};
*/
    /* open/closed configuration
    open-left   open-right
    _            _      ]   [
    _|          |_      [   ]
 part[0] is top, part[1] is bottom
 Target 2 is open-left above open-right, distractor 5 is open-right above open-left
 */
    Shape_parts_t perc_spat_config{spat_configNull};    // open-right/left spatial configuration
    static constexpr Shape_parts_t spat_configNull = {'.','.'};
    static constexpr Shape_parts_t spat_config2 = {'L','R'}; 
    static constexpr Shape_parts_t spat_config5 = {'R','L'};
    static constexpr Shape_parts_t spat_configE = {'R','R'};
    static constexpr Shape_parts_t spat_config3 = {'L','L'};
    static constexpr Shape_parts_t spat_configct = {'R','.'};
    static constexpr Shape_parts_t spat_configcb = {'.','R'};
    static constexpr Shape_parts_t spat_configot = {'L','.'};
    static constexpr Shape_parts_t spat_configob = {'.','L'};

	static std::map<Shape_parts_t, char> parts_to_shape_map;



    // counter for generating unique ID numbers
    static int objectnum;

    // default ctor needed
    Search_object();
    // full constructor sets object name and location as supplied
    Search_object(const char * name, GU::Point location_, GU::Size size_,
        char color, char shape, char orientation, bool is_target = false);  // these are true physical properties
    // partial constructor allows later setting of object ID and location
    Search_object(const char * name, GU::Size size_, char color, char shape, char orientation, bool is_target = false);
    // set object ID
    void set_ID();
    // set object location
    void set_location(GU::Point location_);
    
    void print_sens_spat_config(std::ostream& os) const;
    void print_perc_spat_config(std::ostream& os) const;


};

// used throughout program
using Search_objects_t = std::vector<Search_object> ;

std::ostream& operator<< (std::ostream& os, const Search_object& so);

#endif /* SEARCH_OBJECT_H */


------------------------------------------------------------
// Search_object.cpp
//  
//
//  Created by David Kieras on 3/7/18.
//

#include "Search_object.h"
#include "Program_constants.h"
//#include "EPICLib/Random_utilities.h"
#include "Random_utilities.h"
#include <iostream>

using namespace std;

// counter for object names
int Search_object::objectnum = 0;


std::map<Shape_parts_t, char> Search_object::parts_to_shape_map = {
		{spat_configNull, '.'},
		{spat_config2, '2'},
		{spat_config5, '5'},
		{spat_configE, 'E'},
		{spat_config3, '3'},
		{spat_configct, 'c'},
		{spat_configcb, 'c'},
		{spat_configot, 'o'},
		{spat_configob, 'o'}
		};


// default ctor
Search_object::Search_object() :
    ID(-1)
{
}

// full ctor
Search_object::Search_object(const char * name_, GU::Point location_, GU::Size size_, char color_, char shape_, char orientation_, bool is_target_) :
    name(name_), ID(++objectnum), size(size_), location(location_),
    phys_color(color_), phys_orientation(orientation_), phys_shape(shape_), is_target(is_target_),
    sens_color('.'), sens_orientation('.'), sens_shape('.'),
    perc_color('.'), perc_orientation('.'), perc_shape('.')
    {}
// partial ctor to allow setting ID and location later
Search_object::Search_object(const char * name_, GU::Size size_, char color_, char shape_, char orientation_, bool is_target_) :
    name(name_), size(size_),
    phys_color(color_), phys_orientation(orientation_), phys_shape(shape_), is_target(is_target_),
    sens_color('.'), sens_orientation('.'), sens_shape('.'),
    perc_color('.'), perc_orientation('.'), perc_shape('.')
    {}

void Search_object::set_ID()
{
    ID = ++objectnum;
}

void Search_object::set_location(GU::Point location_)
{
    location = location_;
}

void Search_object::print_sens_spat_config(ostream& os) const
{
    for(int i = 0; i < perc_spat_config.size(); i++)
        os << sens_spat_config[i];
}
void Search_object::print_perc_spat_config(ostream& os) const
{
    for(int i = 0; i < perc_spat_config.size(); i++)
        os << perc_spat_config[i];
}


ostream& operator<< (ostream& os, const Search_object& obj)
{
    os << obj.name << ' ' << obj.ID << ' ' << obj.location << ' ' << obj.eccentricity << ' ' << obj.size << ' '
        << obj.nearest_neighbor_distance << ' ';
        // output crowding status
    os << ((obj.current_crowding) ? 'c' : 'u') << ' '
        << ((obj.known_perceptual_properties) ? 'k' : ' ') << ' ';
    os << obj.phys_color << obj.phys_orientation << obj.phys_shape << " -> "
        << obj.sens_color << obj.sens_orientation << obj.sens_shape << ' ';
    obj.Search_object::print_sens_spat_config(os);
   //     os << " -> " << obj.perc_color << obj.perc_orientation << obj.perc_shape << ' ';
        os << " -> " << ((obj.known_perceptual_properties) ? "[" : "") << obj.perc_color << obj.perc_orientation << obj.perc_shape << ' ';
    obj.Search_object::print_perc_spat_config(os);
    	os << ((obj.known_perceptual_properties) ? "]" : "");
    // below added circa 4/26/19 to go with VM2b
    os << ' ';
    os << obj.Shape_information_gain << ' ';
    os << obj.Feature_information_gain;

    return os;
}



------------------------------------------------------------
// Standard_utility_symbols.h
#ifndef STANDARD_UTILITY_SYMBOLS_H
#define STANDARD_UTILITY_SYMBOLS_H

#include "Symbol.h"

// These are standard names known throughout the architecture as symbols
// in fundamental utility functions and containers.
// The string-representation is normally the same as the name with the "_c" suffix removed.
// See Standard_utility_symbols.cpp for the actual definition.
// As new names are used in the architecture, they should be added to this set,
// and the code should use these variable names instead of explicit strings.

extern const Symbol Default_c;
extern const Symbol Absent_c;
extern const Symbol Unknown_c;
extern const Symbol None_c;
extern const Symbol Nil_c;
extern const Symbol Empty_string_c;

#endif

------------------------------------------------------------
// Standard_utility_symbols.cpp
#include "Standard_utility_symbols.h"

// These are standard names known throughout the architecture as symbols
// in fundamental utility functions and containers.
// The string-representation is normally the same as the name with the "_c" suffix removed.
// See Standard_utility_symbols.cpp for the actual definition.
// As new names are used in the architecture, they should be added to this set,
// and the code should use these variable names instead of explicit strings.

const Symbol Default_c("Default");
const Symbol Absent_c("Absent");
const Symbol Unknown_c("Unknown");
const Symbol None_c("None");
const Symbol Nil_c;
const Symbol Empty_string_c("");


------------------------------------------------------------
// Statistics.h
#ifndef STATISTICS_H
#define STATISTICS_H

#include <vector>
#include <limits>

/*
Use these classes to accumulate running statistics values easily.
reset() - zero the internal variables
update() - add a new data value, updating current average
get_n, mean - return the current values
*/

class Mean_accumulator{
public:
	Mean_accumulator()
		{
			reset();
		}
	void reset()
		{
			n = 0;
			total = 0.0;
			total2 = 0.0;
            min_value = std::numeric_limits<double>::max();
            max_value = std::numeric_limits<double>::min();
		}
	long get_n() const
		{return n;}
	double get_mean() const
		{return (n > 0) ? (total / n) : 0.0;}
	long get_total() const
		{return total;}
	double get_sample_var() const;
	double get_sample_sd() const;
	double get_est_var() const;
	double get_est_sd() const;
	double get_sdm() const;
    double get_half_95_ci() const;
    double get_min() const {return min_value;}
    double get_max() const {return max_value;}
	void update(double x)
		{
			total += x;
			total2 += x*x;
			n++;
            if(x < min_value) min_value = x;
            if(x > max_value) max_value = x;
		}
    void operator() (double x)
        {update(x);}
			
private:
	long n;
	long double total;
	long double total2;
    double min_value;
    double max_value;
};


// This class provides a proportion of the calls to update with a true argument.
// reset() - zero the internal variables
// update() - add a new data value, updating current statistics
// get_count, _n, _proportion - return the current values

class Proportion_accumulator{
public:
	Proportion_accumulator()
		{
			reset();
		}
	void reset()
		{
			count = 0;
			n = 0;
		}
	
	long get_count() const
		{return count;}
	long get_n() const
		{return n;}
	double get_proportion() const
		{return (n > 0) ? double(count) / n : 0.;}
			
	void update(bool count_it)
		{
			n++;
			if(count_it)
				count++;
		}
    void operator() (bool count_it)
        {update(count_it);}
private:
	long count;
	long n;
};

// Accumulate data values into bins and provide the proportion in each bin.
// Initialize with the number of bins and the size in each bin.
// The first bin starts at 0.  Values too big or too small are accumulated
// in the smallest or largest bin.
// Output for( i = 0; i < n_bins; i++) then i*bin_size is the upper bound on values in that bin
// so for bins 0, 25, 50, count for first is number < 0; second is number >= 0 and < 25, third is >= 25 and < 50
class Distribution_accumulator{
public:
	Distribution_accumulator(int n_bins_, double bin_size_, double baseline_ = 0.0) :
        n_bins(n_bins_), bin_size(bin_size_), baseline(baseline_)
		{
			reset();
		}
	void reset()
		{
			n = 0;
            bins.clear();
            bins.resize(n_bins);
            min_value = std::numeric_limits<double>::max();
            max_value = std::numeric_limits<double>::min();
		}

    // These accessor provide individual bins, or the whole distribution
    int get_n_bins() const
        {return n_bins;}
    double get_bin_size() const
        {return bin_size;}
    double get_baseline() const
        {return baseline;}
	long get_n() const
		{return n;}
//    double get_bin_upper_bound(int bin) const
//        {return bin_size * (bin + 1);}
    double get_bin_lower_bound(int bin) const
        {return bin_size * bin + baseline;}
    double get_bin_upper_bound(int bin) const
        {return bin_size * (bin + 1) + baseline;}
	double get_min() const
		{return min_value;}
	double get_max() const
		{return max_value;}
    int get_bin_count(int bin) const
        {
            return bins[bin];
        }
    double get_bin_proportion(int bin) const
        {
            return (n) ? double(bins[bin]) / n : 0.0;
        }
	// subscript operator returns the proportion
    double operator[] (int bin) const
        {
            return (n) ? double(bins[bin]) / n : 0.0;
        }
    std::vector<double> get_distribution() const
        {
            std::vector<double> result(n_bins);
            for(int i = 0; i < n_bins; i++)
                result[i] = (n) ? double(bins[i]) / n : 0.0;
            return result;
        }
    
    void update(double x)
		{
			n++;
/*            if(x < 0.) {
                bins[0]++;
                return;
                }
            int i = int(x / bin_size);
            if(i < 0)
                bins[0]++;
            else if(i >= n_bins)
                bins[n_bins - 1]++;
            else
                bins[i]++;
*/
            int i = int((x - baseline)/ bin_size);
            if(i < 0)
                bins[0]++;
            else if(i >= n_bins)
                bins[n_bins - 1]++;
            else
                bins[i]++;

            if(x < min_value)
                min_value = x;
            if(x > max_value)
                max_value = x;
		}

    void operator() (double x)
        {update(x);}
			
    // add the other's counts to this accumulator, update min/max values
    // (n_bins, bin_size must be the same)
    void add_counts(const Distribution_accumulator& other);
    void update(const Distribution_accumulator& other)
    	{add_counts(other);}
			
private:
    int n_bins;
    double bin_size;
    double baseline;
	long n;
    double min_value;
    double max_value;
    std::vector<int> bins;
};

// Accumulate integer data values into bins and provide the proportion in each bin.
// Each integer value corresponds to one bin.
// Initialize with the number of bins and the first bin value (defaults to 0).
// The number of bins = last bin value - first bin value + 1.
// The value of the last bin = number of bins - first bin value - 1.
// If count_out_of_range_into_bins is true, values too big or too small are accumulated in the smallest or largest bin;
// if false, they are not included in the bin but are still counted as part of the total n for the distribution.
// Output for( i = 0; i < n_bins; i++) access the number or proportion in each bin

class Discrete_distribution_accumulator{
public:
	Discrete_distribution_accumulator(int n_bins_, int first_bin_value_ = 0, bool count_out_of_range_into_bins = false);

    void reset()
		{
			n = 0;
            bins.clear();
            bins.resize(n_bins);
            min_value = std::numeric_limits<int>::max();
            max_value = std::numeric_limits<int>::min();
		}

    // These accessor provide individual bins, or the whole distribution
    int get_n_bins() const
        {return n_bins;}
    int get_bin_size() const
        {return bin_size;}
    int get_first_bin_value() const
        {return first_bin_value;}
    int get_last_bin_value() const
        {return last_bin_value;}
	long get_n() const
		{return n;}
	int get_min() const
		{return min_value;}
	int get_max() const
		{return max_value;}
    int get_bin_count(int bin) const
        {
            return bins[bin];
        }
    double get_bin_proportion(int bin) const
        {
            return (n) ? double(bins[bin]) / n : 0.0;
        }
	// subscript operator returns the proportion
    double operator[] (int bin) const
        {
            return (n) ? double(bins[bin]) / n : 0.0;
        }
    std::vector<double> get_distribution() const
        {
            std::vector<double> result(n_bins);
            for(int i = 0; i < n_bins; i++)
                result[i] = (n) ? double(bins[i]) / n : 0.0;
            return result;
        }
    
    void update(int x)
		{
			n++;
            if(x < min_value)
                min_value = x;
            if(x > max_value)
                max_value = x;
            // normalize x by subtracting first_bin_value
            int i = x - first_bin_value;
            if(i < 0) {
				if(count_out_of_range_into_bins) {
					bins[0]++;
					}
                }
            else if(i >= n_bins) {
				if(count_out_of_range_into_bins) {
					bins[n_bins - 1]++;
					}
				}
            else
                bins[i]++;
		}

    void operator() (double x)
        {update(x);}
			
    // add the other's counts to this accumulator, update min/max values
    // (n_bins, bin_size must be the same)
    void add_counts(const Discrete_distribution_accumulator& other);
    void update(const Discrete_distribution_accumulator& other)
    	{add_counts(other);}
			
private:
    int n_bins;
    int first_bin_value;
    int last_bin_value;
	bool count_out_of_range_into_bins;
    int bin_size = 1; // always one
    long n;
    int min_value;
    int max_value;
    std::vector<int> bins;
};

// This class manages a vector of Proportion_accumulators; get_bin(int bin) returns
// a reference to the Proportion_accumulator bin.
// Given data pairs (x, v), where v is true or false, accumulate the
// proportion of true values for different values of x classifed into bins
// This is different from a distribution because each bin contains Proportion_accumulator,
// not just a simple count. The total number of data pairs is given by the total number of
// true/false cases counted in all of the bins.
class Binned_proportion_accumulators {
public:
	Binned_proportion_accumulators(int n_bins_, double bin_size_) :
        n_bins(n_bins_), bin_size(bin_size_)
		{
			reset();
		}
	void reset()
		{
			n = 0;
            bins.clear();
            bins.resize(n_bins);
		}

    int get_n_bins() const
        {return n_bins;}
    double get_bin_size() const
        {return bin_size;}
	// get the total count in all bins - number of cases supplied
	int get_n() const
		{
			int sum = 0;
			for(const auto& bin : bins) {sum += bin.get_n();}
			return sum;
		}
    // These accessors provide access to individual proportion accumulators
	const Proportion_accumulator& get_bin(int bin) const
		{return bins[bin];}
	const Proportion_accumulator& operator[] (int bin) const
		{return bins[bin];}
	void update(double x, bool v)
		{
			int i = int(x / bin_size);
            if(i < 0)
                bins[0].update(v);
            else if(i >= n_bins)
                bins[n_bins - 1].update(v);
            else
                bins[i].update(v);;
		}
     void operator() (double x, bool v)
        {update(x, v);}
    
private:
    int n_bins;
    double bin_size;
	long n;
    std::vector<Proportion_accumulator> bins;
};

// This class manages a vector of Mean_accumulators; get_bin(int bin) returns
// a reference to the Mean_accumulators bin.
// Given data pairs (x, v), accumulate the mean, etc, for the value v
// for different values of x classifed into bins
// This is different from a distribution because each bin contains Mean_accumulator,
// not just a simple count. The total number of data pairs is given by the total number of
// data values supplied to all of the bins.

class Binned_mean_accumulators {
public:
	Binned_mean_accumulators(int n_bins_, double bin_size_) :
        n_bins(n_bins_), bin_size(bin_size_)
		{
			reset();
		}
	void reset()
		{
            bins.clear();
            bins.resize(n_bins);
		}

    // These accessor provide individual bins, or the whole distribution
    int get_n_bins() const
        {return n_bins;}
    double get_bin_size() const
        {return bin_size;}
    double get_bin_upper_bound(int bin) const
        {return bin_size * bin + 1;}
	// get the total count in all bins - number of cases supplied
	int get_n() const
		{
			int sum = 0;
			for(const auto& bin : bins) {sum += bin.get_n();}
			return sum;
		}
	// get the minimum value over all cells
	double get_min() const
		{
			double min_value = std::numeric_limits<double>::max();
			for(const auto& bin : bins) {
				if(bin.get_min() < min_value) min_value = bin.get_min();
				}
			return min_value;
		}
	// get the maximum value over all cells
	double get_max() const
		{
			double max_value = std::numeric_limits<double>::min();
			for(const auto& bin : bins) {
				if(bin.get_max() > max_value) max_value = bin.get_max();
				}
			return max_value;
		}
	const Mean_accumulator& get_bin(int bin) const
		{return bins[bin];}
	const Mean_accumulator& operator[] (int bin) const
		{return bins[bin];}
	void update(double x, double v)
		{
			int i = int(x / bin_size);
            if(i < 0)
                bins[0].update(v);
            else if(i >= n_bins)
                bins[n_bins - 1].update(v);
            else
                bins[i].update(v);;
		}
     void operator() (double x, double v)
        {update(x, v);}
			
	void update_with_bin_means(const Binned_mean_accumulators& other);
	void update_with_bin_proportions(const Binned_proportion_accumulators& other);
	
private:
    int n_bins;
    double bin_size;
    std::vector<Mean_accumulator> bins;
};

// Accumulate data for a correlation coeficient and regression line
// Like the others, this class uses the one-pass approach which
// can be numerically unreliable under some conditions
class Correl_accumulator {
public:
	Correl_accumulator()
		{reset();}
	void reset()
		{
			n = 0;
			sumx = 0.;
			sumy = 0.;
			sumxy = 0.;
			sumx2 = 0.;
			sumy2 = 0.;
		}

	void update(double x, double y)
		{
			n++;
			sumx += x;
			sumy += y;
			sumxy += x*y;
			sumx2 += x*x;
			sumy2 += y*y;
		}

     void operator() (double x, double y)
        {update(x, y);}
			
	int get_n() const
		{return n;}

	double get_r() const;

	double get_slope() const
		{
			double numerator = n * sumxy - sumx * sumy;
			double denominator = n * sumx2 - sumx * sumx;
			return (denominator > 0.) ? numerator / denominator : 0.0;
		}

	double get_intercept() const
		{
			return (n) ? (sumy - get_slope() * sumx) / n : 0.0;
		}

	double get_rsq() const
		{
			double r = get_r();
			return r*r;
		}

private:
	int n;
	long double sumx;
	long double sumy;
	long double sumxy;
	long double sumx2;
	long double sumy2;
};

// Give this class object a series of predicted and observed values,
// and then get the goodness-of-fit metrics for them
// using regression fit and simple average absolute error
class PredObs_accumulator {
public:
	PredObs_accumulator() : sum_abs_error(0.), sum_error_prop(0.), total_error2(0.)
		{}
	void reset()
		{
			corr.reset();
            sum_abs_error = 0.;
			sum_error_prop = 0.;
			total_error2 = 0.;
		}

	void update(double predicted, double observed);
	void operator() (double predicted, double observed)
        {update(predicted, observed);}
    
	int get_n() const
		{return corr.get_n();}

	double get_rsq() const
		{return corr.get_rsq();}

	double get_slope() const
		{return corr.get_slope();}

	double get_intercept() const
		{return corr.get_intercept();}

    double get_avg_abs_error() const
        {return (sum_abs_error / corr.get_n());}

	double get_avg_abs_rel_error() const
		{return (sum_error_prop / corr.get_n()) * 100.;}

	double get_rmse() const;

private:
	Correl_accumulator corr;
    long double sum_abs_error;
	long double sum_error_prop;
    long double total_error2;
};



#endif

------------------------------------------------------------
// Statistics.cpp
#include "Statistics.h"
#include <cmath>
#define USE_ASSERT_THROW

#ifdef USE_ASSERT_THROW
#include "Assert_throw.h"
#define ASSERT Assert
#endif

#ifdef USE_CASSERT
#include <cassert>
#define ASSERT assert
#endif

using namespace std;


double Mean_accumulator::get_sample_var() const
{
	if(n < 2)
		return 0.0;
	long double mean = get_mean();
	
	return (total2 / n - mean * mean);
}

double Mean_accumulator::get_sample_sd() const
{
	return sqrt(get_sample_var());
}

double Mean_accumulator::get_est_var() const
{
	return (n < 2) ? 0.0 : (double(n) / (n -1)) * get_sample_var();
}

double Mean_accumulator::get_est_sd() const
{
	return sqrt(get_est_var());
}

double Mean_accumulator::get_sdm() const
{
	return (n < 2) ? 0.0 : get_est_sd() / sqrt(n);
}

double Mean_accumulator::get_half_95_ci() const
{
	return (n < 2) ? 0.0 : 1.96 * get_sdm();
}

void Distribution_accumulator::add_counts(const Distribution_accumulator& other)
{
    ASSERT(n_bins == other.n_bins);
    ASSERT(bin_size == other.bin_size);
    for(int i = 0; i < n_bins; i++)
        bins[i] += other.bins[i];
    n += other.n;
	if(other.min_value < min_value)
		min_value = other.min_value;
	if(other.max_value > max_value)
		max_value = other.max_value;
}

Discrete_distribution_accumulator::Discrete_distribution_accumulator(int n_bins_, int first_bin_value_, bool count_out_of_range_into_bins_) :
        n_bins(n_bins_), first_bin_value(first_bin_value_), last_bin_value(n_bins_ - first_bin_value_ - 1),
		count_out_of_range_into_bins(count_out_of_range_into_bins_)
{
    ASSERT(n_bins >= 2);
    ASSERT(last_bin_value >= first_bin_value + 1);
    reset();
}

void Discrete_distribution_accumulator::add_counts(const Discrete_distribution_accumulator& other)
{
	ASSERT(n_bins == other.n_bins);
    ASSERT(first_bin_value == other.first_bin_value);
    ASSERT(last_bin_value == other.last_bin_value);
    ASSERT(count_out_of_range_into_bins == other.count_out_of_range_into_bins);
	for(int i = 0; i < n_bins; i++)
        bins[i] += other.bins[i];
    n += other.n;
	if(other.min_value < min_value)
		min_value = other.min_value;
	if(other.max_value > max_value)
		max_value = other.max_value;
}


void Binned_mean_accumulators::update_with_bin_means(const Binned_mean_accumulators& other)
{
	ASSERT(n_bins == other.get_n_bins() && bin_size == other.get_bin_size());
	for(int i = 0; i < n_bins; i++) {
		bins[i].update(other[i].get_mean());
		}
}

void Binned_mean_accumulators::update_with_bin_proportions(const Binned_proportion_accumulators& other)
{
	ASSERT(n_bins == other.get_n_bins() && bin_size == other.get_bin_size());
	for(int i = 0; i < n_bins; i++) {
		bins[i].update(other[i].get_proportion());
		}
}


double Correl_accumulator::get_r() const
{
	long double numerator = n * sumxy - sumx * sumy;
	long double denominator1 = n * sumx2 - sumx * sumx;
	long double denominator2 = n * sumy2 - sumy * sumy;
	long double r = (denominator1 > 0. && denominator2 > 0.) ?
		numerator / (sqrt(denominator1) * sqrt(denominator2)) : 0.0;
	return r;
}


void PredObs_accumulator::update(double predicted, double observed)
{
    corr.update(predicted, observed);
    long double error = predicted - observed;
    total_error2 += error * error;
    sum_abs_error += fabs(error);
    // compute proportion of absolute error relative to observed value
    sum_error_prop += fabs(error) / observed;
}

double PredObs_accumulator::get_rmse() const
{
    return sqrt(total_error2 / corr.get_n());
}



------------------------------------------------------------
// Strategy9c.cpp
//  SimpleCrowdingModel
//
//  Created by David Kieras on 4/9/18.
//
// revised 7/13/19 to include recording of mean error RT

#include "Strategy.h"
#include "Display_creator.h"
//#include "Visual_modelV2a.h"
#include "Visual_modelV2e.h"
#include "Trial_statistics.h"
//#include "EPICLib/Random_utilities.h"
#include "Random_utilities.h"


#include <iostream>
#include <iomanip>
#include <sstream>
#include <set>
#include <cassert>

using namespace std;

// implements unified parameterized strategy pre-5a, 6b, 7
// version 9c allows search objects to be reordered by Visual_model
string Strategy::get_strategy_id_string()
{
    ostringstream oss;
    oss << "Strategy 9c with";
//    oss << "Strategy 9b look at furthest with";
    if(!confirm_positive_c && !confirm_negative_c)
        oss << " no confirmation";
    else {
        oss << " confirm";
        if(confirm_positive_c && confirm_negative_c)
            oss << " positive & negative";
        else if(confirm_positive_c)
            oss << " positive";
        else if(confirm_negative_c)
            oss << " negative";
        }
    if(max_n_fixations_g == 0)
        oss << ", eyes fixed";
    else if(max_n_fixations_g > 0 && max_n_fixations_g < 99)
 //       oss << ", respond absent if no target present within " << max_n_fixations_g << " fixations";
        oss << ", respond absent if no target present within " << calculate_average_max_n_fixations() << " fixations";
    else
        oss << ", search until target found";
    if(too_many_illusory_targets_c < 99)
        oss << ", or if " << too_many_illusory_targets_c << " illusory targets found";
    oss << "\n";
    if(confirm_positive_c)
        oss << "Skip confirm positive if fixation within " << close_enough_to_bypass_confirm_positive_c << " DVA of apparent target. ";
    if(confirm_negative_c)
    oss << "Skip confirm negative if at least " << enough_fixations_to_bypass_confirm_negative_c << " fixations already made. ";
    oss << "\n";

    return oss.str();
}


// input is time in ms, output is that time rounded up to next cycle start
double time_at_next_cycle_start(double time);


/* implements unified parameterized strategy pre-5a, 6b, 7 */
void Strategy::do_trial()
{
    // flip a coin to decide whether to use max_n_fixations_g + 0 or + 1, but only if max_n_fixations_g < 99
    int trial_max_n_fixations = max_n_fixations_g;
    if(max_n_fixations_mixture_probability_c > 0.0 && max_n_fixations_g < 99 && biased_coin_flip(max_n_fixations_mixture_probability_c))
        trial_max_n_fixations = max_n_fixations_g + 1;

    // bool target_nominated = false; // used to be here, moved to top of search loop 4/26/19
    
    n_repeat_fixations = 0;
    n_target_fixations = 0;
    fixated_objects.clear(); // keeps track of which objects have been fixated

    int n_illusory_target_fixations = 0;
    double stimulus_onset_time = uniform_random_variable(500., 50.); // start trial at a jittered time in range 450 - 550
    double time = stimulus_onset_time; // time for execution from stimulus presentation to response
    double motor_response_time = 0; // time for motor response
    double total_time = 0.;
    
    search_objects = display_creator.create_display(condition, set_size, polarity);

    Visual_model visual_model(condition, search_objects);
    // time allowance for initial registration in sensory store
    time = time + visual_delay_time_g;
    // round up to next rule execution cycle
    time = time_at_next_cycle_start(time);
    time += rule_execution_time_c; // for startup rule execution


    while(true) { // start of search loop
    if(strategy_trace_output) {
        cout << "Fixation: " << visual_model.get_n_fixations() << endl;
        }

    // 2/3/20 Visual_model can reorder the objects for crowding computation
    GU::Point eye_location = visual_model.get_eye_location();
    
    // this sort order is used to order nominated objects, since linear search will find the first match
    // as of 2/3/20 this ordering of nominations is closest object first -- least to greatest eccentricity from current eye location.
    sort(search_objects.begin(), search_objects.end(),
        [eye_location](const auto& obj1, const auto& obj2){return(GU::cartesian_distance(obj1.location, eye_location) < GU::cartesian_distance(obj2.location, eye_location));});
//        [eye_location](const auto& obj1, const auto& obj2){return(GU::cartesian_distance(obj1.location, eye_location) > GU::cartesian_distance(obj2.location, eye_location));});

    if(strategy_trace_output && false) {
        cout << trial_g << " Strategy: objects ordered by eccentricity:\n";
        for(const auto& obj : search_objects)
            cout << obj << endl;
        }
    
//    Search_objects_t::iterator nominated_object_it;
    //  initialize this here to reset on every loop iteration
    bool target_nominated = false;

    
    // NOMINATION PHASE nominate next object to look at; either a target or the closest object with unknown properties
    // all objects nominated are the closest for the type of the nomination because the contained has been sorted by eccentricity
    // The nomination results in an iterator, nominated_object_it, being set to either an apparent target, a possible target, or .end()
    switch(condition) {
        case Condition_e::CSF: {
            // object with target color -- apparent target
            if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return obj.perc_color == 'R';} ))
                    != search_objects.end()) {
                target_nominated = true;
                break;
                }
            // object with unknown color -- a possible target
            else if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                        [](const auto& obj){return obj.perc_color == '.';} ))
                        != search_objects.end()) {
                break;
                }
            else {
                assert(nominated_object_it == search_objects.end());
                break;
                }
            }
        case Condition_e::SHP: {
            // object with target shape - an apparent target
            if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return obj.perc_shape == '2';} ))
                    != search_objects.end()) {
                target_nominated = true;
                break;
                }
            // object with unknown shape or mis-recognized shape -- a possible target
            else if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
 //                   [](const auto& obj){return (obj.perc_shape == '.' || (obj.perc_shape != '2' && obj.perc_shape != '5'));} ))
                    [](const auto& obj){return (obj.perc_shape != '2' && obj.perc_shape != '5');} ))
                    != search_objects.end()) {
                break;
                }
             else {
                assert(nominated_object_it == search_objects.end());
                break;
                }
           }
        case Condition_e::COC: {
            // object with both target color and target orientation -- an apparent target
            if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return (obj.perc_color == 'R') && (obj.perc_orientation == 'V');} ))
                    != search_objects.end()) {
                target_nominated = true;
                break;
                }
            // object with target color but unknown orientation  -- a possible target
            else if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return (obj.perc_color == 'R') && (obj.perc_orientation == '.');} ))
                    != search_objects.end()) {
                 break;
                 }
            // object with unknown color but target orientation  -- a possible target
            else if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return (obj.perc_color == '.') && (obj.perc_orientation == 'V');} ))
                    != search_objects.end()) {
                break;
                }
            // object with both unknown color and unknown orientation  -- a possible target
            else if((nominated_object_it = find_if(search_objects.begin(), search_objects.end(),
                    [](const auto& obj){return (obj.perc_color == '.') && (obj.perc_orientation == '.');} ))
                    != search_objects.end()) {
                break;
                }
            else {
                assert(nominated_object_it == search_objects.end());
                break;
                }
            }
        }
     
    // even if no nominations made, it takes a cycle
    time += rule_execution_time_c; // for nomination rules Nominate_color_unknown_CSF, Nominate_match_color_CSF, etc

    // CHOOSE PHASE - act on the nominations and number of fixations

    // at least one choice rule always fires and takes a cycle
    time += rule_execution_time_c;    // for "lookat" or "response" rule
    
    // DECIDE WHAT TO DO NEXT

    // if no search fixations allowed (eyes fixed strategy) then respond depending on target_nominated and terminate the trial
    if(trial_max_n_fixations == 0) {
        if(target_nominated) {
            motor_response_time = make_positive_response();
            }
        else {
            motor_response_time = make_negative_response();
            }
        break; // out of search loop
        }
    assert(trial_max_n_fixations > 0);

    // if no nominated next_obj, we need to consider a negative response
    if(nominated_object_it == search_objects.end()) {
        // if not confirming negative responses, make the negative response and terminate the trial
        if(!confirm_negative_c) {
            motor_response_time = make_negative_response();
            break; // out of search loop
            }
        // check to see if maximum number of fixations have already been made; if so, respond negative
        if(visual_model.get_n_fixations() >= trial_max_n_fixations) {
            motor_response_time = make_negative_response();
            break; // out of search loop
            }

        // if confirming negative, consider doing a check eye movement
        // restore these special cases 7/13/18
        // go ahead and respond absent if we have made at least one fixation
        if(visual_model.get_n_fixations() >= enough_fixations_to_bypass_confirm_negative_c) {
            motor_response_time = make_negative_response();
            break; // out of search loop
            }

        // do confirmation: move eye to furthest object and check again
        // and either continue the search loop, or make a negative response
        else {
    	    // objects are in order of increasing eccentricity, so last object is furthest away
		    auto& object_to_fixate = search_objects.back();

/*            // objects are in order of decreasing eccentricity, so first object is furthest away
            auto& next_obj = search_objects.front();
*/
//            int saved_object_ID = next_obj.ID;
            if(strategy_trace_output) {
                cout << trial_g << " Strategy: Confirm-Negative look at: " << object_to_fixate << endl;
                }

            double eye_movement_time = move_eye_to_object(visual_model, object_to_fixate);
            time += eye_movement_time;
            // round up to start of next cycle
            time = time_at_next_cycle_start(time);
            time += rule_execution_time_c;    // for eyemovement_done rule execution

            // confirm that it is not a target
            auto& fixated_object = *fixated_object_it;
            bool non_target_confirmed = false;
            switch(condition) {
                case Condition_e::CSF: {
                    if(fixated_object.perc_color == 'G') {
                        non_target_confirmed = true;
                        }
                    break;
                    }
                case Condition_e::SHP: {
                    if(fixated_object.perc_shape == '5') {
                        non_target_confirmed = true;
                        }
                    break;
                    }
                case Condition_e::COC: {
                    if(!((fixated_object.perc_color == 'R') && (fixated_object.perc_orientation == 'V'))) {
                        non_target_confirmed = true;
                        }
                    break;
                    }
                }
            if(non_target_confirmed) {
                motor_response_time = make_negative_response();
                break; // out of search loop - trial is done
                }
            else {
                // NOTE NOTE NOTE ??? if it actually is a target, we waste a nomination cycle!
                // non-target was not confirmed, continue the loop
                continue; // continue search loop
                }
            } // end of case where no fixations made yet        } // end of no-nomination case
        } // end of case where no objects nominated
        
    // here only if an object was nominated:
    assert(nominated_object_it != search_objects.end());
    auto&  nominated_object = *nominated_object_it;
 //   int fixated_object_ID = nominated_object.ID;

    // if a target was not nominated, move the eye, and continue the loop.
    if(!target_nominated) {
        // check to see if maximum number of fixations have already been made; if so, respond negative
        if(visual_model.get_n_fixations() >= trial_max_n_fixations) {
            motor_response_time = make_negative_response();
            break; // out of search loop
            }
        if(strategy_trace_output) {
            cout << trial_g << " Strategy: look at possible target: " << nominated_object << endl;
            }
        double eye_movement_time = move_eye_to_object(visual_model, nominated_object);
        time += eye_movement_time;
        // round up to start of next cycle
        time = time_at_next_cycle_start(time);
        time += rule_execution_time_c;    // for eyemovement_done rule execution
        continue; // continue the search loop
        }
    // if a target was nominated, if confirmation needed, check it out and respond positive or continue the search loop to get a new nomination
    if(target_nominated) {
        // if no confirmation required, go ahead and make the response
        // debug 2/3/20
 //       cout << " *** TARGET NOMINATED: " << nominated_object.name << ' ' << nominated_object.ID << endl;
//        for(const auto& obj : search_objects)
 //           cout << obj << endl;

        if(!confirm_positive_c) {
            motor_response_time = make_positive_response();
            break; // out of search loop - trial is done
            }
        if(nominated_object.eccentricity <= close_enough_to_bypass_confirm_positive_c) {
            // debug 2/3/20
//            cout << " *** CONFIRMATION BYPASSED: " << nominated_object.name << ' ' << nominated_object.ID << endl;
//            for(const auto& obj : search_objects)
//                cout << obj << endl;
            motor_response_time = make_positive_response();
            break; // out of search loop - trial is done
            }
        // else move eye to the nominated object and confirm, but only if more fixations allowed
        // if more fixations are not allowed, respond negative, because target has not in fact been found. 190221
        if(visual_model.get_n_fixations() >= trial_max_n_fixations) {
//            motor_response_time = make_positive_response();  // was in Strategy9
            motor_response_time = make_negative_response();
            break; // out of search loop
            }
        if(strategy_trace_output) {
            cout << trial_g << " Strategy: Confirm-Positive look at apparent target: " << nominated_object << endl;
            }
        // move eye to nominated object
        double eye_movement_time = move_eye_to_object(visual_model, nominated_object);
        time += eye_movement_time;
        // round up to next rule cycle time for confirmation/discomfirmation rule
        time = time_at_next_cycle_start(time);
        time += rule_execution_time_c;    // for confirmation/disconfirmation rule execution
        // fixated_object_it is the valid iterator pointing to the fixated object after processing done
        auto fixated_object = *fixated_object_it;
        // confirm that fixated_object now has target properties
        // debug 2/3/20 NOTE THAT WE DO NOT REORDER THE OBJECTS AT THIS POINT - THEY ARE NOW IN THE ORDER VM2 left them in!
        bool target_confirmed = false;
        switch(condition) {
            case Condition_e::CSF: {
                if(fixated_object.perc_color == 'R') {
                    target_confirmed = true;
                    }
                break;
                }
            case Condition_e::SHP: {
                if(fixated_object.perc_shape == '2') {
                    target_confirmed = true;
                    }
                break;
                }
            case Condition_e::COC: {
                if((fixated_object.perc_color == 'R') && (fixated_object.perc_orientation == 'V')) {
                    target_confirmed = true;
                    }
                break;
                }
            }
        if(target_confirmed) {
            // debug 2/3/20
//            cout << " *** TARGET CONFIRMED: " << fixated_object.name << ' ' << fixated_object.ID << endl;;
//            for(const auto& obj : search_objects)
//                cout << obj << endl;
            motor_response_time = make_positive_response();
            break; // out of search loop - trial is done
            }
        else {
            // target was not confirmed, continue the loop
            // debug 2/3/20
            if(strategy_trace_output) cout << trial_g << " Strategy: apparent target " << fixated_object.name << ' ' << fixated_object.ID << " was not a target" << endl;
//            cout << "Strategy: apparent target " << next_obj.name << ' ' << next_obj.ID << " was not a target" << endl;
            n_illusory_target_fixations++;
            // ??? DK 3/16/19 stop and respond negative if we have seen too many illusory targets
            if(n_illusory_target_fixations >= too_many_illusory_targets_c) {
                motor_response_time = make_negative_response();
                break;  // out of search loop - trial is done
                }
            target_nominated = false; // looks like only place where target_nominated would have to be reset 4/26/19
            continue; // continue search loop
            }
        } // end of target nominated case
    } // end of search loop
    assert(time > 0);
    total_time = (time - stimulus_onset_time) + motor_response_time;
    assert(total_time > 0);

	if(strategy_trace_output || false)
		cout << trial_g << " Strategy: Times: pure visual delay: " << visual_delay_time_g // included in rule time total
		<< ", rules: " << time
		<< ", motor: " << motor_response_time
        << ", total: " << total_time
		<< ", n_fixations: " << visual_model.get_n_fixations()
        << ", n_target_illusory_distractors: " << visual_model.get_n_target_illusory_distractors()
        << ", n_target_illusory_blanks: " << visual_model.get_n_target_illusory_blanks()
		<< ", n_illusory_targets: " << visual_model.get_n_illusory_targets()
//		<< ' ' << ", n_target_sens_property_present: " << visual_model.get_n_target_sens_property_present()
//		<< ' ' << ", n_target_perc_property_present: " << visual_model.get_n_target_perc_property_present()
//        << ", n_target_property_overwritten: " << visual_model.get_n_target_property_overwritten()
		<< endl;

	
    // response has been made, accumulate statistics
    // accumulate statistics
    if(!error) {
        statistics.mean_RT.update(total_time);
        statistics.mean_n_fixations.update(visual_model.get_n_fixations());
        }
	else {
		statistics.mean_error_RT.update(total_time);
        statistics.mean_error_n_fixations.update(visual_model.get_n_fixations());
		}
    // average illusory target regardless of positive vs negative trials ??? DK 3/13/19 ??? correct only would be more useful here - should be impossible in CSF-
    statistics.mean_illusory_target_fixations.update(n_illusory_target_fixations);
    statistics.mean_repeat_fixations.update(n_repeat_fixations);
    statistics.mean_target_fixations.update(n_target_fixations);
    // calculate how many times the display was processed, including intial onset - counts the initial central fixation
    if(strategy_trace_output)
		cout << visual_model.get_n_fixations()
			<< ' ' << visual_model.get_n_target_illusory_distractors()
			<< ' ' << visual_model.get_n_target_illusory_blanks()
			<< ' ' << visual_model.get_n_illusory_targets()
//			<< ' ' << visual_model.get_n_target_sens_property_present()
//			<< ' ' << visual_model.get_n_target_perc_property_present()
//			<< ' ' << visual_model.get_n_target_property_overwritten()
			<< endl;
    double n_process_display = 1 + visual_model.get_n_fixations();
    statistics.mean_n_target_illusory_distractor.update(visual_model.get_n_target_illusory_distractors() / n_process_display);
    statistics.mean_n_target_illusory_blank.update(visual_model.get_n_target_illusory_blanks() / n_process_display);
    statistics.mean_n_illusory_targets.update(visual_model.get_n_illusory_targets() / n_process_display);
//    statistics.mean_n_target_sens_property_present.update(visual_model.get_n_target_sens_property_present() / n_process_display);
//    statistics.mean_n_target_perc_property_present.update(visual_model.get_n_target_perc_property_present() / n_process_display);
//    statistics.mean_n_target_property_overwritten(visual_model.get_n_target_property_overwritten() / n_process_display);
    return;
} // end of do_trial


double Strategy::move_eye_to_object(Visual_model& visual_model, Search_object& object_to_fixate)
{
    int fixated_object_ID = object_to_fixate.ID; // save ID for later, since iterators may become invalid
    // move eye - nominated_object_it is pointing to next object to look at
    double eye_movement_time = visual_model.move_eye(object_to_fixate);
    // Visual_model may reorder the objects, so the iterator is now invalid
    // reset it using the saved ID.
    fixated_object_it = find_if(search_objects.begin(), search_objects.end(),
        [fixated_object_ID] (Search_object& obj) {return obj.ID == fixated_object_ID;});
    assert(fixated_object_it != search_objects.end());
    Search_object& fixated_obj = *fixated_object_it;
    if(fixated_obj.is_target) // it might be the target we are looking at
        n_target_fixations++;
    auto result = fixated_objects.insert(fixated_object_ID);
    if(!result.second) // ID was not inserted because already there
        n_repeat_fixations++;
    return eye_movement_time;
}

// return the response time
double Strategy::make_positive_response()
{
    bool processing_error = (polarity == Polarity_e::NEGATIVE); // error if correct response was opposite to positive
    bool oops_error = biased_coin_flip(positive_baseline_error_rate_c);
    error = processing_error || oops_error; // error rate is at least guess_error rate;
    response = (!error) ? true : false; // if no error, response is yes - true -; opposite if error
    // debug 02/03/20
    if(strategy_trace_output) {
        cout << trial_g << " Strategy: target present";
        if(processing_error)
            cout << " *** FALSE ALARM ERROR";
        else if(oops_error)
            cout << " *** OOPS FALSE ALARM ERROR";
        cout  << endl;
        }
    statistics.prop_errors.update(error);
//        statistics.prop_illusory_target_fixations.update(false); //??? What is baseline for proportion denominator?
    double motor_response_time = (response == last_response) ? same_key_motor_response_time_c : different_key_motor_response_time_c;
    last_response = response;
    return motor_response_time;
}

double Strategy::make_negative_response()
{
    bool processing_error = (polarity == Polarity_e::POSITIVE); // error if correct response was opposite to positive
    bool oops_error = biased_coin_flip(negative_baseline_error_rate_c);
    error = processing_error || oops_error; // error rate is at least guess_error rate;
    response = (!error) ? false : true; // if no error, response is no - false -; opposite if error
    if(strategy_trace_output) {
        cout << trial_g << " Strategy: target absent";
       if(processing_error)
            cout << " *** MISS ERROR";
        else if(oops_error)
            cout << " *** OOPS MISS ERROR";
        cout << endl;
        }
    statistics.prop_errors.update(error);
    double motor_response_time = (response == last_response) ? same_key_motor_response_time_c : different_key_motor_response_time_c;
    last_response = response;
    return motor_response_time;
}

// input is time in ms, output is that time rounded up to next cycle start
double time_at_next_cycle_start(double time)
{
    int integral_part = int(time);
    int n_cycles = integral_part / rule_execution_time_c;
    double result = (n_cycles + 1) * rule_execution_time_c;
    assert(result > 0.);
    return result;
}

------------------------------------------------------------
// Symbol.h

/*
Symbol objects hold a character string value and/or a numeric value.
These values can be replaced, but not modified inside the Symbol.
If a string value is present, it is used as the basis for comparison.
The numeric value is compared only if no string value is present.

The numeric value is either a single double value, a Point, or a vector of Points.
The single value is stored in both the x and y values of a Point.
Numeric values are compared using the vector comparisons regardless of whether they
are a single value, a single Point, or a vector of Points.

If both string and numeric values are present, the Symbol serves as a "named value" - especially
useful if the value is a vector of many Points.

The strings and numeric values are stored as unique reference counted objects. New values
are compared to the stored values, and re-used if they match. 

Creating, copying, assigning, destroying increment/decrement the reference count for
a value; if the reference count goes to zero, the string or number object id deleted.

The whole goal of this class is to enable copy & assignment to be done with shallow copy and 
equality comparison to be done with only pointer comparison. Ordering comparison
is slowed down by the indirection. Since the strings and numeric value cannot be modified, 
no time is wasted on copying values, but the lookup for previous matching values slows down
creation of new values.

When the symbol is initialized or assigned from a non-symbol string, a check can made about
whether the string can be parsed completely as a number (using std::strtod). If it can, then 
the symbol is set up as containing a numeric value only. If not, it is set up as containing
a string. 

Any attempt to extract a value from a Symbol that does not correspond to the stored value
results in a Symbol_exception being thrown. 

Extracting a std::string always succeeds - the resulting string contains either the string
value of the symbol or the numeric value as written with the output operator.
*/

#ifndef SYMBOL_H
#define SYMBOL_H

#include "Point.h"
#include "Exception.h"
#include "Symbol_memory.h"

#include <iosfwd>
#include <string>	// for interface with strings
#include <vector>	// for interface with vectors
#include <list>		// for typedefs and operators for lists of Symbols

class Symbol_exception : public Exception {
public:
	Symbol_exception(const char * in_msg) : Exception(in_msg) {}
};

namespace GU = Geometry_Utilities;

//class Symbol_memory_Str_rep;
//class Symbol_memory_Vec_rep;

class Symbol {
public:
// change to default to check for number always 6/8/03
// it would be better if these constructors were all explicit to encourage 
// defining const Symbols as much as possible for speed, but there are many
// cases were the string is likely to be unique. Explicit remomved 7/30/06
	Symbol(const char * c = "Nil", bool check_for_number = true)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(c, check_for_number);}

	Symbol(const std::string& s, bool check_for_number = true)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(s.c_str(), check_for_number);}

	Symbol(double in_x)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(in_x);}

	Symbol(int in_x)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(double(in_x));}

	Symbol(long in_x)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(double(in_x));}

	Symbol(double in_x, double in_y)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(GU::Point(in_x, in_y));}

	Symbol(GU::Point p)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(p);}

	Symbol(const std::vector<GU::Point>& v)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(v);}

	Symbol(const char * c, const std::vector<GU::Point>& v)
		: str_rep_ptr(0), vec_rep_ptr(0)
		{setup(c, false); setup(v);}

	~Symbol();

	Symbol(const Symbol& src);
	// uses copy-swap
	const Symbol& operator= (const Symbol& rhs);
	// assign from a C-string
	const Symbol& operator= (const char * rhs);
	const Symbol& operator= (double rhs);
	// swap the members
	void swap(Symbol& other);
	
	bool has_string_value() const
		{return str_rep_ptr;}
	bool has_numeric_value() const
		{return vec_rep_ptr;}
	bool has_single_numeric_value() const;
	bool has_Point_numeric_value() const;
	bool has_pair_numeric_value() const;
	bool has_multiple_numeric_value() const;

	// the accessors will throw an exception if expected data is not present
	// string value must be present for these
	char const * c_str() const;
	int size() const;
	int length() const;
	std::string str() const; 
	// some form of numeric value must be present for these
	double get_numeric_value() const;	// must be a single double value
	GU::Point get_Point() const;	// must be a single Point
	// get x or y part of a pair of values
	double get_x() const;	// must be a single Point present
	double get_y() const;	// must be a single Point present
	// get the vector of Points
	const std::vector<GU::Point>& get_Point_vector() const;	// must be any numeric value


/* Comparisons
Symbol1			Symbol2				Comparison
string only		string only			strings are compared
string only		number only			number is always less-than string
string & number	string				strings are compared, number ignored
string & number	string & number		strings are compared, numbers ignored
string & number	number only			number is always less-than string & number
number only		number only			numbers are compared

Number comparisons are done with vector and Point comparison rules

Comparisons with non-Symbol types follow the corresponding rules 
*/

	bool operator== (const Symbol& rhs) const;
	bool operator!= (const Symbol& rhs) const;
	bool operator< (const Symbol& rhs) const;
	bool operator<= (const Symbol& rhs) const; 
	bool operator> (const Symbol& rhs) const ;
	bool operator>= (const Symbol& rhs) const;

	bool operator== (const char * rhs) const;	
	bool operator!= (const char * rhs) const;		
	bool operator<  (const char * rhs) const;	
	bool operator<= (const char * rhs) const;	
	bool operator>  (const char * rhs) const;	
	bool operator>= (const char * rhs) const;	
	bool operator== (double rhs) const;
	bool operator!= (double rhs) const;
	bool operator<  (double rhs) const;
	bool operator<= (double rhs) const;
	bool operator>  (double rhs) const;
	bool operator>= (double rhs) const;
	
	// non-member comparisons
	friend bool operator== (const char * lhs, const Symbol& rhs);
	friend bool operator!= (const char * lhs, const Symbol& rhs);
	friend bool operator<  (const char * lhs, const Symbol& rhs);
	friend bool operator<= (const char * lhs, const Symbol& rhs);
	friend bool operator>  (const char * lhs, const Symbol& rhs);
	friend bool operator>= (const char * lhs, const Symbol& rhs);
	friend bool operator== (double lhs, const Symbol& rhs);
	friend bool operator!= (double lhs, const Symbol& rhs);
	friend bool operator<  (double lhs, const Symbol& rhs);
	friend bool operator<= (double lhs, const Symbol& rhs);
	friend bool operator>  (double lhs, const Symbol& rhs);
	friend bool operator>= (double lhs, const Symbol& rhs);

	friend std::ostream& operator<< (std::ostream&, const Symbol&);
			
private:
	Symbol_memory_Str_rep * str_rep_ptr;	// zero if no string value is present
	Symbol_memory_Vec_rep * vec_rep_ptr;	// zero if no numeric value is present
	//	str_rep_ptr	vec_rep_ptr		contents of this Symbol
	//	!= 0		0				string only - comparisons done on strings
	//	0			!= 0			numeric value(s) only - comparisons done on values
	//	!= 0		!= 0			"named values" - comparisons done on strings
	//	0			0				not possible

	// setup a Symbol from a supplied c-string
	void setup(const char * c, bool check_for_number);
	// check for the string containing a number, set up this symbol if so.
	void setup(double x);
	void setup(GU::Point p);
	void setup(const std::vector<GU::Point>& v);
};

// following inlined here for speed given possible uses in containers and algorithms
inline bool Symbol::operator== (const Symbol& rhs) const
{	
	// to be equal, both pointers must be equal
	// equal if pointing to the same value or are zero	
	return (str_rep_ptr == rhs.str_rep_ptr) && (vec_rep_ptr == rhs.vec_rep_ptr);
}

inline bool Symbol::operator!= (const Symbol& rhs) const
{
	return !(this->operator== (rhs));
}

std::ostream& operator<< (std::ostream&, const Symbol&);

typedef std::list<Symbol> Symbol_list_t;


#endif

------------------------------------------------------------
// Symbol.cpp

// This defines a class of read-only strings suitable for use as symbols.
// Once created, they can be only output, compared, or destroyed. 
// They can be initialized with or compared to C strings.
// A const pointer to the internal C string can be obtained.

// When the Symbol object is constructed from a string that can be parsed as a number,
// the numeric value is saved and a flag is set the the object contains a numeric value.
// If the object is constructed from a number, a default string is created for it using
// default iostream output rules.
// The numeric value can be accessed, and used in comparisons, but may not be changed.

#include "Symbol.h"
#include "Symbol_memory.h"
#include "Assert_throw.h"
#include "Utility_templates.h"

#include <iostream>
#include <sstream>
#include <cstring>
#include <cstdio>
#include <utility>	// for swap

using std::ostream; using std::cerr; using std::endl;
using std::set;
using std::ostringstream;
using std::string;
using std::strlen; using std::strcpy; using std::strtod; using std::strcmp;
using std::size_t;
using std::vector;

//#include "Output_tee.h"

// setup an Symbol from a supplied c-string
// check for it being a number, save it as a number if so
// if not a number, enter the c-string;

// needs fixing to allow for named value!!!

void Symbol::setup(const char * c, bool check_string_for_number)
{
	int len = int(strlen(c));	// calculate this only once
	// if len is zero, we must want an empty string
	// and it can't contain anything else
	if(len) {
		if (check_string_for_number) {
			// check for the c-string containing a number
			// the whole string must be the number
			char * ptr;
			double temp = strtod(c, &ptr);
			if (ptr == c + len) {
				str_rep_ptr = 0;	// we have a number
				setup(temp); // set it up instead
				return;		// setup is done
				}
			}
		}
	// if check not asked for, this is not a number, store as string only
	// set up the string representation in memory
	str_rep_ptr = Symbol_memory::get_instance().find_or_insert(c, len);
}


void Symbol::setup(double x_)
{
	vector<GU::Point> v(1, GU::Point(x_, x_));
	// set up a single value
	vec_rep_ptr = Symbol_memory::get_instance().find_or_insert(true, v);
}
	
void Symbol::setup(GU::Point p)
{
	vector<GU::Point> v(1, p);
	vec_rep_ptr = Symbol_memory::get_instance().find_or_insert(false, v);
}
	
void Symbol::setup(const vector<GU::Point>& v_)
{
	vec_rep_ptr = Symbol_memory::get_instance().find_or_insert(false, v_);
}

		
Symbol::~Symbol()
{
	// if this symbol has a string, remove it if it is the last
	if(str_rep_ptr) {
		Symbol_memory::get_instance().remove_if_last(str_rep_ptr);
		}
	// ditto if it has a vector
	if(vec_rep_ptr) {
		Symbol_memory::get_instance().remove_if_last(vec_rep_ptr);
		}
}

Symbol::Symbol(const Symbol& src)
{
	str_rep_ptr = src.str_rep_ptr;
	if(str_rep_ptr)
		(str_rep_ptr->count)++;
	vec_rep_ptr = src.vec_rep_ptr;
	if(vec_rep_ptr)
		(vec_rep_ptr->count)++;
}
	
// uses copy-swap
const Symbol& Symbol::operator= (const Symbol& rhs) 
{
	Symbol temp(rhs);
	swap(temp);
	return *this;
}
		
// swap the members; note how reference counts of reps are unchanged
void Symbol::swap(Symbol& other)
{
	std::swap(str_rep_ptr, other.str_rep_ptr);
	std::swap(vec_rep_ptr, other.vec_rep_ptr);
}

// assign from a C-string
const Symbol& Symbol::operator= (const char * rhs)
{
	Symbol temp(rhs);
	swap(temp);
	return *this;
}

const Symbol& Symbol::operator= (double rhs)
{
	Symbol temp(rhs);
	swap(temp);
	return *this;
}

bool Symbol::has_single_numeric_value() const 
{
	return vec_rep_ptr && vec_rep_ptr->single;
}

bool Symbol::has_Point_numeric_value() const 
{
	return vec_rep_ptr && !(vec_rep_ptr->single) && (vec_rep_ptr->vec.size() == 1);
}

bool Symbol::has_pair_numeric_value() const 
{
	return has_Point_numeric_value();
}

bool Symbol::has_multiple_numeric_value() const 
{
	return vec_rep_ptr && (vec_rep_ptr->vec.size() > 2);
}


// return a C string or throw an exception if there is none
char const * Symbol::c_str() const
{
	if (str_rep_ptr)
		return str_rep_ptr->cstr;
	throw Symbol_exception("No C-string for Symbol");
}	

int Symbol::size() const
{
	if (str_rep_ptr)
//		return int(strlen(str_rep_ptr->cstr));
		return int(str_rep_ptr->cstr_len);
	throw Symbol_exception("No C-string size for Symbol");
}

int Symbol::length() const
{
	if (str_rep_ptr)
		return int(str_rep_ptr->cstr_len);
	throw Symbol_exception("No C-string length for Symbol");
}

// return a string containing either the C string or a string
// representation of the numeric value using the output operator
string Symbol::str() const
{
	ostringstream oss;
		oss << *this;
	return oss.str();	
}

// get the numeric values or throw an exception if not present
double Symbol::get_numeric_value() const
{
	if(has_single_numeric_value())
		return vec_rep_ptr->vec[0].x;
	throw Symbol_exception("Symbol is not single numeric value");
}

GU::Point Symbol::get_Point() const
{
	if(has_Point_numeric_value())
		return vec_rep_ptr->vec[0];
	throw Symbol_exception("Symbol is not a Point numeric value");
}

double Symbol::get_x() const
{
	if(has_pair_numeric_value())
		return vec_rep_ptr->vec[0].x;
	throw Symbol_exception("Symbol is not a pair numeric value");
}

double Symbol::get_y() const
{
	if(has_pair_numeric_value())
		return vec_rep_ptr->vec[0].y;
	throw Symbol_exception("Symbol is not a pair numeric value");
}

const vector<GU::Point>& Symbol::get_Point_vector() const
{
	if(has_numeric_value())
		return vec_rep_ptr->vec;
	throw Symbol_exception("Symbol is not a vector value");
}
	
/* Comparisons
Symbol1			Symbol2				Comparison
string only		string only			strings are compared
string only		number only			number is always less-than string
string & number	string				strings are compared, number ignored
string & number	string & number		strings are compared, numbers ignored
string & number	number only			number is always less-than string & number
number only		number only			numbers are compared

Number comparisons are done with vector and Point comparison rules

Comparisons with non-Symbol types follow the corresponding rules 
*/

// not inlined because Str_rep and Vec_rep would have to be visible
bool Symbol::operator< (const Symbol& rhs) const 
{
	// if the two are different string configurations, can quit immediately
	if(str_rep_ptr && !rhs.str_rep_ptr) {
		//Assert(rhs.vec_rep_ptr);
		return false;	// str is not less than number only
		}
	else if(!str_rep_ptr && rhs.str_rep_ptr) {
		//Assert(vec_rep_ptr);
		return true;	// number only is less than str
		}
	else if(str_rep_ptr && rhs.str_rep_ptr)	{ // if both strings, compare them
		if(str_rep_ptr == rhs.str_rep_ptr) {	// same strings 
			return false;	// same str is not less than
			}
		// compare the different strings for less than
		else 
			return (strcmp(str_rep_ptr->cstr, rhs.str_rep_ptr->cstr) < 0);
		}
	// no strings, just numbers
	//Assert(vec_rep_ptr && rhs.vec_rep_ptr);
	return vec_rep_ptr->vec < rhs.vec_rep_ptr->vec;
}


// define these in terms of operator<
bool Symbol::operator<= (const Symbol& rhs) const 
{
	return  !(rhs < *this);
}

bool Symbol::operator> (const Symbol& rhs) const 
{
	return rhs < *this;
}
		
bool Symbol::operator>= (const Symbol& rhs) const 
{
	return  !(*this < rhs);
}

/* Member comparisons */
bool Symbol::operator== (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) == 0) : false;
}

bool Symbol::operator!= (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) != 0) : true;
}

bool Symbol::operator<  (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) < 0) : true;
}

bool Symbol::operator<=  (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) <= 0) : true;
}

bool Symbol::operator>  (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) > 0) : false;
}

bool Symbol::operator>=  (const char * rhs) const
{
	return (str_rep_ptr) ? (std::strcmp(str_rep_ptr->cstr, rhs) >= 0) : false;
}

bool Symbol::operator== (double rhs) const
{
	return (str_rep_ptr) ? false : (get_numeric_value() == rhs);
}

bool Symbol::operator!= (double rhs) const
{
	return (str_rep_ptr) ? true : (get_numeric_value() != rhs);
}

bool Symbol::operator<  (double rhs) const
{
	return (str_rep_ptr) ? false : (get_numeric_value() < rhs);
}

bool Symbol::operator<=  (double rhs) const
{
	return (str_rep_ptr) ? false : (get_numeric_value() <= rhs);
}

bool Symbol::operator>  (double rhs) const
{
	return (str_rep_ptr) ? true : (get_numeric_value() > rhs);
}

bool Symbol::operator>=  (double rhs) const
{
	return (str_rep_ptr) ? true : (get_numeric_value() >= rhs);
}

/* Non-member comparisons */

// if rhs has no string, result is false
bool operator== (const char * lhs, const Symbol& rhs)  
//	{return (rhs.str_rep_ptr) ? (strcmp(lhs, rhs.str_rep_ptr->cstr) == 0) : false;}
{
	return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) == 0) : false;}

// if rhs has no string, result is true
bool operator!= (const char * lhs, const Symbol& rhs)  
	{return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) != 0) : true;}

// if rhs has no string, result is false - numbers are less
bool operator< (const char * lhs, const Symbol& rhs)  
	{return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) < 0) : false;}

bool operator<= (const char * lhs, const Symbol& rhs)  
	{return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) <= 0) : false;}

// if rhs has no string, result is true - numbers are less
bool operator> (const char * lhs, const Symbol& rhs)  
	{return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) > 0) : true;}

bool operator>= (const char * lhs, const Symbol& rhs)  
	{return (rhs.has_string_value()) ? (strcmp(lhs, rhs.c_str()) >= 0) : true;}

// if rhs is not a number, result is false
bool operator== (double lhs, const Symbol& rhs)
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() == lhs) : false;}

// if rhs is not a number, result is true
bool operator!= (double lhs, const Symbol& rhs)
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() != lhs) : true;}
		
// if rhs is not a number, result is true - numbers are less
bool operator< (double lhs, const Symbol& rhs) 
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() < lhs) : true;}

bool operator<= (double lhs, const Symbol& rhs) 
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() <= lhs) : true;}

// if rhs is not a number, result is false - numbers are less
bool operator> (double lhs, const Symbol& rhs) 
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() > lhs) : false;}

bool operator>= (double lhs, const Symbol& rhs) 
	{return (rhs.has_numeric_value()) ? (rhs.get_numeric_value() >= lhs) : false;}
	
// output operators
ostream& operator<< (ostream& os, const Symbol& s) 
{
	if (s.has_string_value()) {
		os << s.c_str();
		}
	else if (s.has_single_numeric_value()) {
		os << s.get_numeric_value();
		}
	else if (s.has_pair_numeric_value()) {
		os << '(' << s.get_x() << ", " << s.get_y() << ')';
		}
	// nothing for other cases
	return os;
}

------------------------------------------------------------
// Symbol_memory.h
#ifndef SYMBOL_MEMORY_H
#define SYMBOL_MEMORY_H

#include "Point.h"
#include <vector>
#include <cstring>
#include <set>

namespace GU = Geometry_Utilities;

// A class for keeping reference-counted strings
class Symbol_memory_Str_rep {
public:
	Symbol_memory_Str_rep(long count_, const char * cstr_, int cstr_len_) :
		count(count_), cstr(cstr_), cstr_len(cstr_len_)
		{}
	friend class Symbol;
	friend class Symbol_memory;
	friend struct Symbol_memory_less_Str_rep;
private:
	long count;
	const char * cstr;
	size_t cstr_len;
};

// function object class for ordering Str_reps
struct Symbol_memory_less_Str_rep {
	bool operator() (const Symbol_memory_Str_rep * lhs, const Symbol_memory_Str_rep * rhs) const
		{
			return (std::strcmp(lhs->cstr, rhs->cstr) < 0);
		}
};

// A class for keeping reference-counted vectors of Points
class Symbol_memory_Vec_rep {
public:
	Symbol_memory_Vec_rep(long count_, bool single_, const std::vector<GU::Point>& v_) :
		count(count_), single(single_), vec(v_)
		{}
	friend class Symbol;
	friend class Symbol_memory;
	friend struct Symbol_memory_less_Vec_rep;
private:
	long count;
	bool single;	// true if only a single value (not a pair) is stored
	const std::vector<GU::Point> vec;
};

// function object class for ordering Vec_reps
// single value is less than non-single
struct Symbol_memory_less_Vec_rep {
	bool operator() (const Symbol_memory_Vec_rep * lhs, const Symbol_memory_Vec_rep * rhs) const
		{
			if(lhs->single == rhs->single)
				return (lhs->vec < rhs->vec);
			else if (lhs->single && !(rhs->single))
				return true;
			else
				return false;
		}
};


// This singleton class encapsulates a container for Str_reps and Vec_reps
// It is used by the Symbol class, hence most of the members are private
class Symbol_memory {
public:
	static Symbol_memory& get_instance();

	~Symbol_memory()
		{clear();}

	friend class Symbol;

private:
	static Symbol_memory * Symbol_memory_ptr;
	// privatize the members for this singleton class
	Symbol_memory()
		{}
	Symbol_memory(const Symbol_memory&);
	Symbol_memory& operator= (const Symbol_memory&);
	
	// deallocate all of the memory.
	void clear();
		
	// return the pointer to the Symbol_memory_Str_rep if the cstring is already present, 
	// add the Symbol_memory_Str_rep if it isn't, and return the pointer
	// len is supplied because it is already computed - to save time
	Symbol_memory_Str_rep * find_or_insert(const char * p, int len);
	
	// decrement reference count for pointed-to string, and deallocate it if no
	// longer in use
	void remove_if_last(Symbol_memory_Str_rep *);

	typedef std::set<Symbol_memory_Str_rep *, Symbol_memory_less_Str_rep> Str_rep_ptr_set_t;
	Str_rep_ptr_set_t str_rep_ptr_set;
	std::size_t total_Str_rep_size;
	
	// Return the pointer to the Symbol_memory_Vec_rep if it is already present, 
	// add the Symbol_memory_Vec_rep if it isn't, and return the pointer
	Symbol_memory_Vec_rep * find_or_insert(bool single_, const std::vector<GU::Point>& v_);
			
	// decrement reference count for pointed-to string, and deallocate it if no
	// longer in use
	void remove_if_last(Symbol_memory_Vec_rep *);

	// Not clear that there  is any advantage to keeping a uniquified container of Vec_reps
	// that are shared between Symbols ...
	typedef std::set<Symbol_memory_Vec_rep *, Symbol_memory_less_Vec_rep> Vec_rep_ptr_set_t;
	Vec_rep_ptr_set_t vec_rep_ptr_set;
};

#endif


------------------------------------------------------------
// Symbol_memory.cpp
#include "Symbol_memory.h"
#include "Utility_templates.h"
#include "Assert_throw.h"
#include <cstddef>
#include <algorithm>

using std::vector;
using std::set;
using std::strcpy;
using std::size_t;

Symbol_memory * Symbol_memory::Symbol_memory_ptr = 0;

// A Meyers singleton would get automatically destructed, which is not
// a good idea because the individual Symbols might get destructed afterwards
// everything should get freed except for the final empty container
Symbol_memory& Symbol_memory::get_instance()
{
	if(!Symbol_memory_ptr)
		Symbol_memory_ptr = new Symbol_memory;
	return *Symbol_memory_ptr;
}


// deallocate the memory for the c-strings and the vectors
void Symbol_memory::clear()
{
	for_each(vec_rep_ptr_set.begin(), vec_rep_ptr_set.end(), Delete());
	vec_rep_ptr_set.clear();
	
	for(Str_rep_ptr_set_t::iterator it = str_rep_ptr_set.begin(); it != str_rep_ptr_set.end(); it++) {
		Symbol_memory_Str_rep * p = *it;
		delete[] p->cstr;
		delete p;
		}
	// deallocate the set contents
	str_rep_ptr_set.clear();
}

// return the pointer to the Symbol_memory_Str_rep if the cstring is already present, 
// add the Symbol_memory_Str_rep if it isn't, and return the pointer
// len is supplied because it is already computed - to save time
Symbol_memory_Str_rep * Symbol_memory::find_or_insert(const char * p, int len)
{
	Symbol_memory_Str_rep sr(0, p, 0);
	Str_rep_ptr_set_t::iterator it = str_rep_ptr_set.find(&sr);
	if (it == str_rep_ptr_set.end()) {
		char * cp = new char[len + 1];
		strcpy(cp, p);
		Symbol_memory_Str_rep * str_rep_ptr = new Symbol_memory_Str_rep(1, cp, len);
		str_rep_ptr_set.insert(str_rep_ptr);
		return str_rep_ptr;
		}
	else {
		Symbol_memory_Str_rep * str_rep_ptr = *it;
		// increment the count
		(str_rep_ptr->count)++;
		return str_rep_ptr;
		}
}

Symbol_memory_Vec_rep * Symbol_memory::find_or_insert(bool single_, const std::vector<GU::Point>& v_)
{
	Symbol_memory_Vec_rep vr(0, single_, v_);
	Vec_rep_ptr_set_t::iterator it = vec_rep_ptr_set.find(&vr);
	if (it == vec_rep_ptr_set.end()) {
		Symbol_memory_Vec_rep * vec_rep_ptr = new Symbol_memory_Vec_rep(1, single_, v_);
		vec_rep_ptr_set.insert(vec_rep_ptr);
		return vec_rep_ptr;
		}
	else {
		Symbol_memory_Vec_rep * vec_rep_ptr = *it;
		// increment the count
		(vec_rep_ptr->count)++;
		return vec_rep_ptr;
		}
}

// decrement the reference count; remove it and free the memory if it was the last use
void Symbol_memory::remove_if_last(Symbol_memory_Str_rep * str_rep_ptr)
{
	(str_rep_ptr->count)--;
	if(str_rep_ptr->count == 0) {
		// remove it from the set
		Str_rep_ptr_set_t::iterator it = str_rep_ptr_set.find(str_rep_ptr);
		Assert(it != str_rep_ptr_set.end());
		str_rep_ptr_set.erase(it);
		delete str_rep_ptr;
		}
}

// decrement the reference count; remove it and free the memory if it was the last use
void Symbol_memory::remove_if_last(Symbol_memory_Vec_rep * vec_rep_ptr)
{
	(vec_rep_ptr->count)--;
	if(vec_rep_ptr->count == 0) {
		// remove it from the set
		Vec_rep_ptr_set_t::iterator it = vec_rep_ptr_set.find(vec_rep_ptr);
		Assert(it != vec_rep_ptr_set.end());
		vec_rep_ptr_set.erase(it);
		delete vec_rep_ptr;
		}
}


------------------------------------------------------------
// Symbol_utilities.h
#ifndef SYMBOL_UTILITIES_H
#define SYMBOL_UTILITIES_H

#include "Symbol.h"
#include <iosfwd>

class Output_tee;

// Return the Symbol that is nth in the list, starting with 0
// if not found, or not legal n, an empty Symbol is returned
Symbol get_nth_Symbol(const Symbol_list_t& in_list, int n);

void print_Symbol_list (const Symbol_list_t& in_list);
void print_Symbol_list (const Symbol_list_t& in_list, Output_tee& ot);

Symbol_list_t cstr_to_Symbol_list(const char *);

Symbol int_to_Symbol(int i);

// return a Symbol consisting of the supplied string(s) or Symbol(s) 
// followed by the digits of the supplied integer

Symbol concatenate_to_Symbol(const char * str, long i);
Symbol concatenate_to_Symbol(const Symbol& sym, long i);
Symbol concatenate_to_Symbol(const char * str1, const char * str2, long i);
Symbol concatenate_to_Symbol(const Symbol& sym1, const Symbol& sym2, long i);
Symbol concatenate_to_Symbol(const char * str, const Symbol& sym, long i);
Symbol concatenate_to_Symbol(const Symbol& sym, const char * str, long i);
Symbol concatenate_to_Symbol(const char * str1, const char * str2, const char * str3, long i);
Symbol concatenate_to_Symbol(const char * str1, const Symbol& sym, const char * str3, long i);

// defined in Symbol.h
//typedef std::list<Symbol> Symbol_list_t;

typedef std::list<Symbol_list_t> list_Symbol_list_t;

typedef std::list<list_Symbol_list_t> list_list_Symbol_list_t;

std::string concatenate_to_string(const Symbol_list_t&);

std::ostream& operator<< (std::ostream& os, const Symbol_list_t& symbol_list);
std::ostream& operator<< (std::ostream& os, const list_Symbol_list_t& list_symbol_list);
std::ostream& operator<< (std::ostream& os, const list_list_Symbol_list_t& list_list_symbol_list);

#endif

------------------------------------------------------------
// Symbol_utilities.cpp
#include "Symbol_utilities.h"
#include "Output_tee.h"

#include <iostream>
#include <sstream>
#include <string>


using std::ostream;	using std::cout;	using std::endl;
using std::ostringstream;	using std::istringstream;
using std::string;


// Return the Symbol that is nth in the list, starting with 0
// if not found, or not legal n, an empty Symbol is returned
Symbol get_nth_Symbol(const Symbol_list_t & in_list, int n)
{
	if (n > int((in_list.size() - 1)))
		return Symbol();
	
	Symbol_list_t ::const_iterator it = in_list.begin();
	
	for (int i = 0; i < n; i++)
		it++;
		
	if (it != in_list.end())
		return *it;
	else
		return Symbol();
}

void print_Symbol_list (const Symbol_list_t & in_list)
{
	Symbol_list_t ::const_iterator it;
	
	for (it = in_list.begin(); it != in_list.end(); it++)
		cout << *it << endl;
}

void print_Symbol_list(const Symbol_list_t & in_list, Output_tee& ot)
{
	Symbol_list_t ::const_iterator it;
	
	for (it = in_list.begin(); it != in_list.end(); it++)
		ot << *it << endl;
}

// return a list of Symbols, where each sysbol is the whitespace delimited sequence
// in the input C-string.  e.g. "A B CD E" => (A B CD E)
Symbol_list_t cstr_to_Symbol_list(const char * str)
{
	istringstream is(str);
	string s;
	Symbol_list_t result;
	while(is >> s)
		result.push_back(Symbol(s));
	return result;
}


Symbol int_to_Symbol(int i)
{
	ostringstream ss;
	ss << i;
	return Symbol(ss.str());
}

Symbol concatenate_to_Symbol(const char * str, long i)
{
	ostringstream ss;
	ss << str;
	ss << i;
	return Symbol(ss.str());
}

Symbol concatenate_to_Symbol(const Symbol& sym, long i)
{
	return concatenate_to_Symbol(sym.c_str(), i);
}

Symbol concatenate_to_Symbol(const char * str1, const char * str2, long i)
{
	ostringstream ss;
	ss << str1;
	ss << str2;
	ss << i;
	return Symbol(ss.str());
}

Symbol concatenate_to_Symbol(const Symbol& sym1, const Symbol& sym2, long i)
{
	return concatenate_to_Symbol(sym1.c_str(), sym2.c_str(), i);
}

Symbol concatenate_to_Symbol(const char * str, const Symbol& sym, long i)
{
	return concatenate_to_Symbol(str, sym.c_str(), i);
}

Symbol concatenate_to_Symbol(const Symbol& sym, const char * str, long i)
{
	return concatenate_to_Symbol(sym.c_str(), str, i);
}

Symbol concatenate_to_Symbol(const char * str1, const char * str2, const char * str3, long i)
{
	ostringstream ss;
	ss << str1;
	ss << str2;
	ss << str3;
	ss << i;
	return Symbol(ss.str());
}

Symbol concatenate_to_Symbol(const char * str1, const Symbol& sym, const char * str3, long i)
{
	return concatenate_to_Symbol(str1, sym.c_str(), str3, i);
}

// copy the list of Symbols into a space-delimited string
string concatenate_to_string(const Symbol_list_t& symbol_list)
{
	if (symbol_list.empty())
		return string();
		
	ostringstream ss;
	for(Symbol_list_t::const_iterator it = symbol_list.begin(); it != symbol_list.end(); it++) {
		if (it != symbol_list.begin())
			ss << ' ';
		ss << (*it);
		}
	return ss.str();
}

ostream& operator<< (ostream& os, const Symbol_list_t& symbol_list)
{
	if (symbol_list.empty())
		return os;
		
	os << '(';
	for(Symbol_list_t::const_iterator it = symbol_list.begin(); it != symbol_list.end(); it++) {
		if (it != symbol_list.begin())
			os << ' ';
		os << (*it);
		}
	os << ')';

	return os;
}

ostream& operator<< (ostream& os, const list_Symbol_list_t& list_symbol_list)
{
	if (list_symbol_list.empty())
		return os;
		
	for(list_Symbol_list_t::const_iterator it = list_symbol_list.begin(); it != list_symbol_list.end(); it++) {
		if(it != list_symbol_list.begin())
			os << endl;
		os << (*it);
		}
	return os;
}


ostream& operator<< (ostream& os, const list_list_Symbol_list_t& list_list_symbol_list)
{
	if (list_list_symbol_list.empty())
		return os;
		
	for(list_list_Symbol_list_t::const_iterator it = list_list_symbol_list.begin(); it != list_list_symbol_list.end(); it++) {
		if(it != list_list_symbol_list.begin())
			os << endl;
		os << (*it);
		}
	return os;
}

------------------------------------------------------------
// Trial_statistics.h
//
//  Created by David Kieras on 4/9/18.
//

#ifndef TRIAL_STATISTICS_H
#define TRIAL_STATISTICS_H
//#include "EPICLib/Statistics.h"
#include "Statistics.h"

struct Trial_statistics {
    Mean_accumulator mean_RT;
    Mean_accumulator mean_n_fixations;
    Mean_accumulator mean_error_RT;
	Mean_accumulator mean_error_n_fixations;
    Proportion_accumulator prop_errors;
    Mean_accumulator mean_illusory_target_fixations; // mean number of fixations per trial that are illusory target fixations
    Mean_accumulator mean_repeat_fixations; // mean number of fixations per trial that are repeat fixations
    Mean_accumulator mean_target_fixations; // mean number of fixations per trial that are target fixations
    // These accumulators count the mean (over trials) of the number of objects across fixations that have center properties
    // mean number of objects that are illusory distractors per trial that have distractor perceptual property
    Mean_accumulator mean_n_target_illusory_distractor;
    // mean number of objects that are illusory blanks because target sensory property available
    Mean_accumulator mean_n_target_illusory_blank;
    // mean number of objects that are illusory targets per trial
    Mean_accumulator mean_n_illusory_targets;
    // mean number of objects that have the target sensory property
//    Mean_accumulator mean_n_target_sens_property_present;
    // mean number of objects that have the target sensory property
//    Mean_accumulator mean_n_target_perc_property_present;
    // mean number of fixations in which target property was available but got overwritten (sens but not perc)
//    Mean_accumulator mean_n_target_property_overwritten;
};

#endif /* STATISTICS_H */

------------------------------------------------------------
// Trial_statistics.cpp
//
//  Created by David Kieras on 4/9/18.
//

#include "Trial_statistics.h"

------------------------------------------------------------
// Visual_modelV2e.h
//
//  Created by David Kieras on 3/13/18.
//

#ifndef VISUAL_MODEL_H
#define VISUAL_MODEL_H

//#include "EPICLib/Geometry.h"
#include "Geometry.h"
#include "Program_constants.h"
#include "Search_object.h"
#include <iosfwd>
#include <vector>
#include <array>
#include <string>

namespace GU = Geometry_Utilities;
using Search_objects_t = std::vector<Search_object> ;

// A Visual_model object is created for each trial; statistics are collected for that trial only
class Visual_model {
public:
    Visual_model(Condition_e cond_, Search_objects_t& search_objects_)
        : condition(cond_), search_objects(search_objects_)
        {initialize();}
    // the eye starts at the fixation point; call this for each subsequent eye move
    // returns time in ms for eye movement motor activity to be done
    double move_eye(Search_object& saccade_target_obj);
    
    GU::Point get_eye_location() const
        {return eye_location;}
	
	int get_n_fixations() const
		{return n_fixations;}
        
    // how many times a distractor has target perceptual properties
    int get_n_illusory_targets() const
        {return n_illusory_targets;}

    // how many times a target has distractor properties
    int get_n_target_illusory_distractors() const
        {return n_target_illusory_distractors;}

    // how many times a target has blank property when sensory known
    int get_n_target_illusory_blanks() const
        {return n_target_illusory_blanks;}
/*
    // how many times a target has distractor properties
    int get_n_target_sens_property_present() const
        {return n_target_sens_property_present;}

    // how many times a target has distractor properties
    int get_n_target_perc_property_present() const
        {return n_target_perc_property_present;}

    // how many times no object with target property is available (no actual or illusory target)
    int get_n_target_property_overwritten() const
        {return n_target_property_overwritten;}
*/
    static std::string get_id_string();

	
private:
    Search_objects_t& search_objects; // top level creates this container, referenced here
    Condition_e condition;
    GU::Point eye_location;
    // assumes Visual_model created for each trial, so these counters do not aggregate across trials
	int n_fixations = 0;
    // these counters are incremented by process_search_objects, initially and after each fixation
    int n_illusory_targets = 0; // counts how many times a distractor has target perceptual properties
    int n_target_illusory_distractors = 0; // counts how many times a target has distractor properties
	int n_target_illusory_blanks = 0; // counts how many times a target has a blank property when sensory property known
//	int n_target_sens_property_present = 0; // counts how many times an object has the target sensory property
//	int n_target_perc_property_present = 0;// counts how many times an object has the target perceptual property
//    int n_target_property_overwritten = 0;  // counts how many times target preceptual properties were available for a non-target object but the target does not have the target property
    
    // initialize the eye position and process the first set of objects
    void initialize();
    void process_search_objects();
    double make_saccade(GU::Point new_location);
    void update_sensory_properties();
    void find_crowders();
    void update_perceptual_properties();
    void update_Shape_information_gain();
    void update_Feature_information_gain();
    void scramble_perceptual_properties();
    void update_crowding_status();
    void scramble_perc_properties_by_permutation(Search_object& obj, char Search_object::*perc_prop_ptr);
    void scramble_perc_spat_configs_by_permutation(Search_object& obj);
    void recognize_shape(Search_object& obj, char Search_object::*from_ptr, Shape_parts_t Search_object::*to_ptr);
    void print_objects() const;

    void check_target_status(Condition_e condition, const Search_objects_t& search_objects);
    void check_CSF_target_status(const Search_objects_t& search_objects);
    void check_COC_target_status(const Search_objects_t& search_objects);
    void check_SHP_target_status(const Search_objects_t& search_objects);
};







#endif /* Visual_model_h */

------------------------------------------------------------
// Visual_modelV2e.cpp
//  Created by David Kieras on 5/3/19.
//

/*
Version V2e. This version includes internal options for the order in which objects are processed
for crowding computation. This corrects an odd situation in which objects were processed in the
last order that the Strategy had left the Search_objects container, which doesn't make much sense,
becaue the eye was in a different location.

Version V2d. This version eliminates the Search_object current_crpwding member variable which
was being incorrectly set during computing the crowding group to show that an flanking object was crowded
when it not necessarily was crowded. The size of an object's crowding group > 0 means that it has crowders.
This means the current_crowded member variable is redundant.

Version V2c. This version adds options to enable/disable the features added since the original V2.
The option controls are in the Program_constants module.

Version V2b. This version adds an information-gained score for each object
that gives the expected number of objects whose properties would be revealed if this object is fixated;
the expected number is determined by the availability of the surrounding objects assuming that the
point of fixation is a directly on this object.

Version V2a. This version modifies treatment of objects that are in perceptual store state:
The property can be overwritten by a non-blank property, but not by a blank property.
This means that the scrambling can replace a property in perceptual store state,
but only by an available property - a previously available property is still "sticky"
if object falls into peripheral vision, but become subject to crowding anew. 190219

Version V2. This version attempts to fix the problem that V1 disrupted the memory of
the properties for previously fixated objects.

Version that reorders the objects so that crowding computations are performed in a well-defined order -
previous version performed the computation in whatever order the strategy last left the objects in, which
is semi-haphazard.

NOTE: This module is supposed to be independent of the Strategy model -- the Strategy model can ask for the eye to be moved
to an object's location, but this module computations are based only on the current location of the eye, not on the object
whose location was used to fixate.  So if this module reorders the objects in its container (a vector, for speed), the Strategy
cannot assume that an interator or reference to an object in the vector is still valid - it must re-find the object.
*/

#include "Program_constants.h"
#include "Visual_modelV2e.h"
#include "Search_object.h"
//#include "EPICLib/Random_utilities.h"
#include "Random_utilities.h"
#include <iostream>
#include <sstream>
#include <cassert>
#include <cmath>


using namespace std;

// This should be moved to EPICLib as being generally useful
double normal_cdf(double x, double mean, double sd);

bool LinearConstSD_available(double eccentricity, GU::Size size, double a, double b, double sd);
double LinearConstSD_detection_probability(double eccentricity, GU::Size size, double a, double b, double sd);

enum class Crowding_comp_order_e {Last_order, From_furthest, From_closest, Random};

Crowding_comp_order_e crowding_comp_order = Crowding_comp_order_e::From_closest;

/* This model implements crowding for each object in a way that maximizes crowding
and statistically independently of the number of crowding objects.
*/
string Visual_model::get_id_string()
{
    ostringstream oss;
//   oss << "Visual model V2d";
//    oss << "Visual model V2dx";
//    oss << "Visual model V2dy";
//    oss << "Visual model V2dz";
    oss << "Visual model V2e";
	if(!apply_crowding_effect)
    	oss << " with no crowding";
	else
		switch (crowding_comp_order) {
			case Crowding_comp_order_e::Last_order: {
				oss << " with crowding probability cp if (n_crowders > 0), crowding applied to each object in last order, asymetrically between objects differing in eccentricity.\n"; // V2d original
				break;
				}
			case Crowding_comp_order_e::From_furthest: {
				oss << " with crowding probability cp if (n_crowders > 0), crowding applied to each object from greatest eccentricity, asymetrically between objects differing in eccentricity.\n"; // V2dx
				break;
				}
			case Crowding_comp_order_e::From_closest: {
				oss << " with crowding probability cp if (n_crowders > 0), crowding applied to each object from least eccentricity, asymetrically between objects differing in eccentricity.\n"; // V2dy
				break;
				}
			case Crowding_comp_order_e::Random: {
				oss << " with crowding probability cp if (n_crowders > 0), crowding applied to each object in random order, asymetrically between objects differing in eccentricity.\n";// V2dz
				break;
				}
			default:
				assert(!"undefined Crowding_comp_order");
			};

		oss << "Once object properties available with no crowders, ";
        switch (perceptual_store_scramble_policy) {
            case Perceptual_store_scramble_policy_e::NO_OVERWRITES: {
                 oss << "its properties are no longer scrambled.\n";
               break;}
            case Perceptual_store_scramble_policy_e::NON_BLANK_OVERWRITES: {
                oss << "scrambling can replace with only non-blank property.\n";
                break;}
            case Perceptual_store_scramble_policy_e::ALL_OVERWRITES: {
                oss << "scrambling can still change its properties.\n";
                break;}
            }

    if(unitary_shape_c)
        oss << "Unitary shape.";
    else
        oss << "Shape with " << n_shape_parts_c << " parts.";
    
    if(saccade_noise_present_c)
        oss << " Saccade noise present.";
    else
        oss << " No saccade noise.";

    return oss.str();
}
// initialize the eye position and process the first set of objects
void Visual_model::initialize()
{
    n_fixations = 0;
    eye_location = GU::Point(0., 0.);
    // special case set of eccentricity
    for(auto& obj : search_objects)
    	obj.eccentricity = cartesian_distance(eye_location, obj.location);
	process_search_objects(); // present display before update
	if(processed_visual_trace_output || detailed_visual_trace_output) {
		cout << "*** after update for initial display; eye at: " << eye_location << endl;
		print_objects();
		}
}

// move eye to the new location using noisy saccades, then process the objects
double Visual_model::move_eye(Search_object& obj)
{
//    if(strategy_trace_output)
//        cout << "Visual: before saccade to " << obj.ID << " at " << obj.location << endl;
    double time = make_saccade(obj.location); // make the saccade
    n_fixations++;
    process_search_objects();
//    if(strategy_trace_output)
//        cout << "Visual: after saccade to " << obj.ID << " eccentricity is " << obj.eccentricity << endl;
	if(processed_visual_trace_output || detailed_visual_trace_output) {
		cout << "Visual: after update for eye fixation " << n_fixations << " on object " << obj.ID << "; eye at: " << eye_location << endl;
    	print_objects();
    	}
    return time;
}

double Visual_model::make_saccade(GU::Point new_location)
{
    double gain = .95;
    double s = 0.10;
    double saccade_angle_sd = 1;

    if(!saccade_noise_present_c) {
        gain = 1.0;
        s = 0.0;
        saccade_angle_sd = 0;
        }

    double time = 0;
    
    GU::Polar_vector move_vector(eye_location, new_location);
    // no move if eye is actually already there, essentially
    if(move_vector.r < 0.1){
        move_vector.r = 0.0;
        time = em_time_intercept_c; // ???
        return time;
        }
    assert(move_vector.r > 0.);
    // modify r using the saccade distance parameters - gain and s, sample actual distance
    double base_distance = move_vector.r * gain;
    double sd = base_distance * s;
    double actual_distance = base_distance;
    // if shortfall_distance is 50, .1 times that gives sd of 5, 3 sds out is 15 degrees, saccade of 75 degrees
    // require that actual distance is > 0.
    if(sd > 0.0) {
        do {actual_distance = trimmed_normal_random_variable(base_distance, sd, 3.0);
            } while (actual_distance <= 0.0);
        }

    move_vector.r = actual_distance;

    // modify theta using noise parameter independent of r
    double saccade_angle_SD = GU::to_radians(saccade_angle_sd);
    double actual_angle = move_vector.theta;
    if(saccade_angle_SD > 0.0) {
        actual_angle = trimmed_normal_random_variable(actual_angle, saccade_angle_SD, 3.0);
        }
    move_vector.theta = actual_angle;

    // translate the eye location using the computed move_vector
    eye_location = eye_location + move_vector;
    time = em_time_intercept_c + em_time_slope_c * move_vector.r;
    return time;
}


void Visual_model::process_search_objects()
{
    if(detailed_visual_trace_output) {
        cout << "\nVisual: before update of sensory properties; eye at: " << eye_location << endl;
        print_objects();
        }
    update_sensory_properties();
    if(detailed_visual_trace_output) {
		cout << "Visual: after update of sensory properties; eye at: " << eye_location << endl;
		print_objects();
		}
    update_perceptual_properties();
    if(detailed_visual_trace_output) {
        cout << "Visual: after update of perceptual properties; eye at: " << eye_location << endl;
        print_objects();
        }
    if(apply_crowding_effect) {
		switch (crowding_comp_order) {
			case Crowding_comp_order_e::Last_order: {
				// V2d original
				break;
				}
			case Crowding_comp_order_e::From_furthest: {
				// V2dx sort objects from largest eccentricity to smallest before crowding
				sort(search_objects.begin(), search_objects.end(),
					[this](const auto& obj1, const auto& obj2)
					{return (GU::cartesian_distance(obj1.location, eye_location) > 	GU::cartesian_distance(obj2.location, eye_location));});
				break;
			case Crowding_comp_order_e::From_closest: {
				// V2dy sort objects from smallest eccentricity to largest before crowding
				sort(search_objects.begin(), search_objects.end(),
					[this](const auto& obj1, const auto& obj2)
					{return(GU::cartesian_distance(obj1.location, eye_location) < 	GU::cartesian_distance(obj2.location, eye_location));});
				break;
				}
				}
			case Crowding_comp_order_e::Random: {
				// V2dz put objects in a random order before crowding
				shuffle(search_objects.begin(), search_objects.end(), get_Random_engine());
				break;
				}
			default:
				assert(!"undefined Crowding_comp_order");
			};

        if(detailed_visual_trace_output || crowding_visual_trace_output) {
            cout << "Visual: after reordering of search objects relative to " << eye_location << endl;
            print_objects();
            }
        find_crowders();
        scramble_perceptual_properties();
        if(detailed_visual_trace_output || crowding_visual_trace_output) {
            cout << "Visual: after scramble of perceptual properties; eye at: " << eye_location << endl;
            print_objects();
            }
        // now if all of a currently uncrowded object's properties are known, mark it as known
        update_crowding_status();
        }
    
    // always check tha target status for instrumentation purposes
    check_target_status(condition, search_objects);
    
    if(compute_Shape_information_gain_g) {
        update_Shape_information_gain();
        if(detailed_visual_trace_output) {
            cout << "Visual: after update_Shape_information_gain; eye at: " << eye_location << endl;
            print_objects();
            }
        }
    
    if(compute_Feature_information_gain_g) {
        update_Feature_information_gain();
        if(detailed_visual_trace_output) {
            cout << "Visual: after update_Feature_information_gain; eye at: " << eye_location << endl;
            print_objects();
            }
        }

//    if(processed_visual_trace_output || detailed_visual_trace_output || target_status_trace_output) {
//        check_target_status(condition, search_objects);
//        }
}

void Visual_model::print_objects() const
{
	for(auto& obj : search_objects)
		cout << obj << endl;
}


// set the sensory properties - depend only on property, eccentricity, and size - not crowding
void Visual_model::update_sensory_properties()
{
//    cout << "spat_config2: "; segs_out(spat_config2); cout << "\n";
//    cout << "spat_config5: "; segs_out(spat_config5); cout << "\n";

    for(auto& obj : search_objects) {
        obj.eccentricity = cartesian_distance(obj.location, eye_location);
        obj.sens_color = LinearConstSD_available(obj.eccentricity, obj.size, common_intercept_c, color_slope_g, common_sd_c) ? obj.phys_color : '.';
        obj.sens_orientation = LinearConstSD_available(obj.eccentricity, obj.size, common_intercept_c, orientation_slope_g, common_sd_c) ? obj.phys_orientation : '.';
        if(unitary_shape_c) {
            obj.sens_shape = LinearConstSD_available(obj.eccentricity, obj.size, common_intercept_c, shape_slope_g, common_sd_c) ? obj.phys_shape : '.';
            }
        else {
            // for each segment, sample availability, leave as blank if not present in phys
            // synthesize phys_spat_config here instead of having it constructed in the Search object
            Shape_parts_t phys_spat_config;
            if(obj.phys_shape == '2')
                phys_spat_config = Search_object::spat_config2;
            else if(obj.phys_shape == '5')
                phys_spat_config = Search_object::spat_config5;
            else
                phys_spat_config = Search_object::spat_configNull; // no 2 or 5 shape specified, have null spat_config pattern

            // now sample availability for each segment in the physical spat config
            for(int i = 0; i < obj.sens_spat_config.size(); i++) {
                // using shape_slope_g for each shape_part
                obj.sens_spat_config[i] = (LinearConstSD_available(obj.eccentricity, shape_part_size_c, common_intercept_c, shape_slope_g, common_sd_c)) ?
                    phys_spat_config[i] :
                    '.';
                }
            }
        }
}



// simply move known sensory property to perceptual property - if sensory property is unknown, leave perceptual property alone
void Visual_model::update_perceptual_properties()
{
    // update the object properties for each object
    for(auto& obj : search_objects) {
        // simple idea: we always update perceptual properties according to now available sensory properties
        // with the perceptual properties that are present being sticky.
        // cases:
        // sens perc -> perc
        //  .   .       . (unchanged)
        //  .   x       x (unchanged) sticky property
        //  x   .       x (update now available property)
        //  x   x       x (unchanged)
        //  x   y       x (correct bogus previous property)
        // non-null sens_color unconditionally becomes current perceptual color
        if(obj.sens_color != '.') {
            obj.perc_color = obj.sens_color;
            }
        if(obj.sens_orientation != '.') {
            obj.perc_orientation = obj.sens_orientation;
           }
        if(unitary_shape_c) {
            if(obj.sens_shape != '.') {
                obj.perc_shape = obj.sens_shape;
                }
            }
        else {
            for(int i = 0; i < obj.sens_spat_config.size(); i++) {
                if(obj.sens_spat_config[i] != '.')
                    obj.perc_spat_config[i] = obj.sens_spat_config[i];
                }
        
            // set the perc_shape unitary property based on the perc_spat_config pattern
            // object, from-field, to-field
            recognize_shape(obj, &Search_object::perc_shape, &Search_object::perc_spat_config);
            }
        // other cases are automatic
    }
}



// assuming update_sensory_properties already called to update eccentricity
// go through all search objects and find out which other objects are crowding them : the crowding group
// note that the current object is a member of the crowding group, so size of group is >= 1
// TODO: Fix that - no reason that uncrowded object represented that way - see later code checking for > 1 for crowding group size
// if no crowders, mark the object is_uncrowded unless - sticky property
/*
Normally, if an object is crowded by another object, or crowds another object, we scramble properties with the other object.
If it neither is crowded, nor crowds, any other objects, is become uncrowded and thereafter, its assigned properties will be sticky
and not changed even if it crowds in the future, unless the properties are different while it is currently uncrowded.
???  TODO UPDATE THIS COMMENT
*/
void Visual_model::find_crowders()
{
    // first, set the current_crowding to false, we'll update this while finding all the crowding groups
//    for (auto& obj : search_objects)
//       obj.current_crowding = false;

    // now, for each object, find out if it has any crowders, and mark both this and the crowding objects as being currently crowded
    for(int i = 0; i < search_objects.size(); i++) {
        Search_object& obj1 = search_objects[i];
        obj1.crowder_group.clear();
       double critical_spacing = bouma_fraction_c * obj1.eccentricity; // bouma function
        for(int j = 0; j < search_objects.size(); j++) {
            if(i == j)
                continue;
            Search_object& obj2 = search_objects[j];
            double obj_spacing = cartesian_distance(obj1.location, obj2.location);
            if(obj_spacing <= critical_spacing) {
                // add obj2's index to obj1's crowding group
                obj1.crowder_group.push_back(j);
                // mark that obj1 is currently crowded by another object
//                obj1.current_crowding = true;
                // mark that obj2 is currently crowded by another object
//                obj2.current_crowding = true;
                }
            }
        // the crowder group is empty if this object has no crowders, > 1 if it does
        // put this object's index into the crowding group
        if(obj1.crowder_group.size() > 0) {
            obj1.crowder_group.push_back(i); // put index of this object into the list
            }
       }
}


// assumes that update_perceptual_properties and find_crowders has already been called
// scramble if there are crowders
void Visual_model::scramble_perceptual_properties()
{
    // scramble the object properties for each object
    for(auto& obj : search_objects) {
    
        if(detailed_visual_trace_output) {
            cout << "Visual: " << obj.ID << " has " << obj.crowder_group.size() << " objects in its crowder group {";
            for(int crowder_index : obj.crowder_group) {
                cout << ' ' << search_objects[crowder_index].ID;
                }
            cout << " }" << endl;
            }
        // either no crowders for the object, or this object and at least one other object
        assert(obj.crowder_group.size() == 0 || obj.crowder_group.size() >= 2);
        if(obj.crowder_group.size() > 1) {
            switch(condition) {
                case Condition_e::CSF:
                    if(biased_coin_flip(color_crowding_probability_g))
                        scramble_perc_properties_by_permutation(obj, &Search_object::perc_color);
                    break;
                case Condition_e::COC:
                    if(biased_coin_flip(color_crowding_probability_g))
                        scramble_perc_properties_by_permutation(obj, &Search_object::perc_color);
                    if(biased_coin_flip(orientation_crowding_probability_g))
                        scramble_perc_properties_by_permutation(obj, &Search_object::perc_orientation);
                    break;
                case Condition_e::SHP:
                    if(unitary_shape_c) {
                        if(biased_coin_flip(shape_crowding_probability_g))
                            scramble_perc_properties_by_permutation(obj, &Search_object::perc_shape);
                        }
                    else {
                        if(biased_coin_flip(shape_crowding_probability_g))
                            scramble_perc_spat_configs_by_permutation(obj);
                        }
                    break;
                default:
                    assert(!"unknown condition in update_perceptual_properties");
                    break;
        }

            } // end of crowder_group.size() > 1 case
        } // end of each-object loop
}

// scramble the property pointed to by perc_prop_ptr
// but if either obj, or an object in the crowding group is marked as known_perceptual_properties, do not modify its properties
// but its properties are still scrambled with those of other objects and so might be assigned to some other object
void Visual_model::scramble_perc_properties_by_permutation(Search_object& obj, char Search_object::*perc_prop_ptr)
{
    // crowder_group has indicies for crowding objects
    // current object is included in the list
    // get the properties of objects in the crowding grouping and scramble their order
    vector<char> scrambled_properties;
    for(int i = 0; i < obj.crowder_group.size(); i++) {
        int obj_idx = obj.crowder_group[i];
        if(detailed_visual_trace_output) cout << search_objects[obj_idx].ID << ' ' << 	search_objects[obj_idx].*perc_prop_ptr << endl;
        scrambled_properties.push_back(search_objects[obj_idx].*perc_prop_ptr);
        }
    if(detailed_visual_trace_output) {
        for(auto c : scrambled_properties)
            cout << c << ' ';
        cout << endl;
        }
    shuffle(scrambled_properties.begin(), scrambled_properties.end(), get_Random_engine());
    assert(obj.crowder_group.size() == scrambled_properties.size());
    if(detailed_visual_trace_output) {
        for(auto c : scrambled_properties)
            cout << c << ' ';
        cout << endl;
        }

    // go through the crowder list and assign the property in the scrambled_properties list
    // but do not assign to a blank property to an object with known_perceptual_properties 190219
    for(int i = 0; i < obj.crowder_group.size(); i++) {
        int obj_idx = obj.crowder_group[i];
        Search_object& obj = search_objects[obj_idx];
        // if object is not in perceptual store state, always scramble the properties
        if(!obj.known_perceptual_properties) {
            obj.*perc_prop_ptr = scrambled_properties[i];
            if(detailed_visual_trace_output || (obj.is_target && target_status_trace_output))
                cout << obj.ID << "<-" << obj.*perc_prop_ptr << endl;
                }
        else {  // if the object is in perceptual store state, scramble depending on the policy
            switch (perceptual_store_scramble_policy) {
                case Perceptual_store_scramble_policy_e::NO_OVERWRITES: {
                    // do nothing
                    break;}
               
            case Perceptual_store_scramble_policy_e::NON_BLANK_OVERWRITES: {
                // scramble only if the property is non-blank
                // check: assign if either not known_perceptual_properties OR non-blank property
                if(scrambled_properties[i] != '.') {
                        obj.*perc_prop_ptr = scrambled_properties[i];
                    if(detailed_visual_trace_output || (obj.is_target && target_status_trace_output))
                            cout << obj.ID << "<-" << obj.*perc_prop_ptr << endl;
                        }
                    break;}
                
            case Perceptual_store_scramble_policy_e::ALL_OVERWRITES: {
                // scramble the property regardless
                obj.*perc_prop_ptr = scrambled_properties[i];
                if(detailed_visual_trace_output || (obj.is_target && target_status_trace_output))
                    cout << obj.ID << "<-" << obj.*perc_prop_ptr << endl;
                break;}
					
			default: {
				assert(!"Unknown Perceptual_store_scramble_policy in Visual_model");
				break;}
                } // end of switch
            } // end of else
        } // end of for
}


// scramble the property pointed to by perc_prop_ptr
void Visual_model::scramble_perc_spat_configs_by_permutation(Search_object& obj)
{
    // crowder_group has indicies for crowding objects
    // current object is included in the list
    // get the properties of objects in the crowding grouping and scramble their order
    if(detailed_visual_trace_output) {
        cout << obj.ID << " has " << obj.crowder_group.size() << " objects in its crowder group {";
        for(int crowder_index : obj.crowder_group) {
            cout << ' ' << search_objects[crowder_index].ID;
            }
        cout << " }" << endl;
        }

    // the segment properties are all lumped together because it is assumed that each property has a location as a subpart of the object.
    vector<char> scrambled_properties;
    for(int i = 0; i < obj.crowder_group.size(); i++) {
        int obj_idx = obj.crowder_group[i];
        for(int j = 0; j < n_shape_parts_c; j++) {
             scrambled_properties.push_back(search_objects[obj_idx].perc_spat_config[j]);
             }
//        if(detailed_visual_trace_output) cout << search_objects[obj_idx].ID << ' ' << search_objects[obj_idx].*perc_prop_ptr << endl;
        }
    if(detailed_visual_trace_output) {
        for(auto c : scrambled_properties)
            cout << c << ' ';
        cout << endl;
        }
    shuffle(scrambled_properties.begin(), scrambled_properties.end(), get_Random_engine());
    assert(scrambled_properties.size() == obj.crowder_group.size() * n_shape_parts_c);
    if(detailed_visual_trace_output) {
        for(auto c : scrambled_properties)
            cout << c << ' ';
        cout << endl;
        }

    // go through the crowder list and assign the property in the scrambled_properties list
    // but do not assign blank property to an object with known_perceptual_properties
    // !!! NOTE THIS HAS NOT BEEN UPDATED AS ABOVE 190219
    for(int i = 0; i < obj.crowder_group.size(); i++) {
        int obj_idx = obj.crowder_group[i];
        Search_object& obj = search_objects[obj_idx];
        if(!obj.known_perceptual_properties) {
            if(detailed_visual_trace_output || (obj.is_target && target_status_trace_output)) {
                cout << obj.ID << ' '; obj.print_perc_spat_config(cout);
                cout << " -> ";
                }
            for(int j = 0; j < n_shape_parts_c; j++) {
                int prop_idx = i*n_shape_parts_c + j;
                assert(prop_idx < scrambled_properties.size());
                obj.perc_spat_config[j]= scrambled_properties[prop_idx];
                }
            if(detailed_visual_trace_output || (obj.is_target && target_status_trace_output)) {
                cout << obj.ID << ' '; obj.print_perc_spat_config(cout);
                cout << "\n";
                }
            }
        }
    // update each object's shape property based on segments - do here to simplify matching rules
    // note performed on all objects
    for(int i = 0; i < obj.crowder_group.size(); i++) {
        int obj_idx = obj.crowder_group[i];
        recognize_shape(search_objects[obj_idx], &Search_object::perc_shape, &Search_object::perc_spat_config);
        }
}

// set the shape unitary property based on the spat_config pattern
void Visual_model::recognize_shape(Search_object& obj, char Search_object::*to_ptr, Shape_parts_t Search_object::*from_ptr)
{
     obj.*to_ptr = Search_object::parts_to_shape_map[obj.*from_ptr];
/*
        if(obj.*from_ptr == Search_object::spat_config2)
            obj.*to_ptr = '2';
        else if(obj.*from_ptr == Search_object::spat_config5)
            obj.*to_ptr = '5';
        else if(obj.*from_ptr == Search_object::spat_configE)
            obj.*to_ptr = 'E';
        else if(obj.*from_ptr == Search_object::spat_config3)
            obj.*to_ptr = '3';
        else
            obj.*to_ptr = 'X';
*/
}

// for each object, if its relevant properties are all known, and it is currently uncrowded, mark it as known_perceptual_properties
// the size() of the crowder_group unambiguously shows whether the object is crowded or not.
void Visual_model::update_crowding_status()
{
    for(auto& obj : search_objects) {
//    	if(obj.current_crowding) assert(obj.crowder_group.size() > 1);
//    	if(obj.current_crowding && obj.crowder_group.size() == 0) cout <<  obj.crowder_group.size() << ' ';
//        if(!obj.current_crowding) {
        if(obj.crowder_group.size() == 0) {
//        	assert(obj.crowder_group.size() == 0);
            switch(condition) {
                case Condition_e::CSF:
                    if(obj.perc_color != '.') {
                        obj.known_perceptual_properties = true;
                        }
                    break;
                case Condition_e::COC:
                    if(obj.perc_color != '.' && obj.perc_orientation != '.')
                        obj.known_perceptual_properties = true;
                    break;
                case Condition_e::SHP:
                    if(unitary_shape_c) {
                        if(obj.perc_shape != '.') {
                            obj.known_perceptual_properties = true;
                            }
                        }
                    else {
                        if(all_of(obj.perc_spat_config.begin(), obj.perc_spat_config.end(), [](const char c) {return c != '.';})) {
                            obj.known_perceptual_properties = true;
                            }
                        }
                    break;
                default:
                    assert(!"unknown condition in update_crowding_status");
                    break;
                }
            if(detailed_visual_trace_output && obj.known_perceptual_properties) {
                cout << "Visual: " << obj.ID << " now has known perceptual properties " << endl;
                }
            } // end of currently uncrowded object conditional
        } // end of each object loop
}


// update_Shape_information_gain
// For each object obj:
// Calculate the expected number of objects whose blank shape will become known
// if obj is fixated -- included in this calculation is obj as well as the surrounding
// objects. If an object's Shape property is already known (i.e. perc_shape != '.'),
// then it is not included in the expected value, so this is the "gain" in information
// that is expected by fixating obj. Since obj is expected to have eccentricity = 0
// after the move, if its shape is currently known, it will contribute zero, but if
// currently unknown, will contribute approximately 1.0 to the information_gain score
void Visual_model::update_Shape_information_gain()
{
    assert(unitary_shape_c);    // this function not implemented for non-unitary shape property

    // zero the scores first
    for(auto& obj : search_objects) {
        obj.Shape_information_gain = 0.0;
        }

    for(auto& obj : search_objects) {
        // credit this object obj with information gain if perceptual Shape is unknown
        if(obj.perc_shape == '.')
            // Expected value that Shape will be available assuming 0.0 eccentricity if fixated
            obj.Shape_information_gain += LinearConstSD_detection_probability(0.0, obj.size, common_intercept_c, shape_slope_g, common_sd_c);
        // for every other object whose shape is unknown, determine expected value that its Shape will be available if obj is fixated
        for(auto& other_obj : search_objects) {
            if(other_obj.ID == obj.ID)
                continue;
            if(other_obj.perc_shape != '.') // skip if shape is known
                continue;
            double eccentricity = cartesian_distance(obj.location, other_obj.location);
            obj.Shape_information_gain += LinearConstSD_detection_probability(eccentricity, other_obj.size, common_intercept_c, shape_slope_g, common_sd_c);
            }
        }
}


// update_Feature_information_gain
// For each object obj:
// Calculate the expected number of objects whose blank features will become known
// if obj is fixated -- included in this calculation is obj as well as the surrounding
// objects. The score is calculated over both color and shape so that e.g. twice as much
// credit is given if both features would become available instead of just one.
// If an object's color property is already known (i.e. perc_color != '.'),
// then it is not included in the expected value, so this is the "gain" in information
// that is expected by fixating obj. Since obj is expected to have eccentricity = 0
// after the move, if its color is currently known, it will contribute zero, but if
// currently unknown, will contribute approximately 1.0 to the information_gain score
// Likewise, the orientation propety is computed for each object as well, and added to the information_gain score
void Visual_model::update_Feature_information_gain()
{
    assert(unitary_shape_c);    // this function not implemented for non-unitary shape property

    // zero the scores first
    for(auto& obj : search_objects) {
        obj.Feature_information_gain = 0.0;
        }

    for(auto& obj : search_objects) {
        // credit this object obj with information gain if perceptual Color is unknown
        if(obj.perc_color == '.') {
            // Expected value that Color will be available assuming 0.0 eccentricity if fixated
            obj.Feature_information_gain += LinearConstSD_detection_probability(0.0, obj.size, common_intercept_c, color_slope_g, common_sd_c);
            }
        // ditto for Orientation
        if(obj.perc_orientation == '.') {
            // Expected value that Orientation will be available assuming 0.0 eccentricity if fixated
            obj.Feature_information_gain += LinearConstSD_detection_probability(0.0, obj.size, common_intercept_c, orientation_slope_g, common_sd_c);
            }
        // for every other object whose Color is unknown, determine expected value that its Color will be available if obj is fixated
        // and ditto for Orientation
        for(auto& other_obj : search_objects) {
            if(other_obj.ID == obj.ID) // skip if same object
                continue;
            // we only need to compute eccentricity if we are going to need it, tweak later if useful
            if(other_obj.perc_color == '.') {// if color is unknown
                double eccentricity = cartesian_distance(obj.location, other_obj.location);
                obj.Feature_information_gain += LinearConstSD_detection_probability(eccentricity, other_obj.size, common_intercept_c, color_slope_g, common_sd_c);
                }
             if(other_obj.perc_orientation == '.') {// if orientation is unknown
                double eccentricity = cartesian_distance(obj.location, other_obj.location);
                obj.Feature_information_gain += LinearConstSD_detection_probability(eccentricity, other_obj.size, common_intercept_c, orientation_slope_g, common_sd_c);
                }
           }
        }
}



void Visual_model::check_target_status(Condition_e condition, const Search_objects_t& search_objects)
{
    switch(condition) {
        case Condition_e::CSF:
            check_CSF_target_status(search_objects);
            break;
        case Condition_e::COC:
            check_COC_target_status(search_objects);
            break;
        case Condition_e::SHP:
            check_SHP_target_status(search_objects);
            break;
        default:
            assert(!"unknown condition in check_target_status");
            break;
        }
}
//target_status_trace_output
void Visual_model::check_CSF_target_status(const Search_objects_t& search_objects)
{
//    bool target_sens_property_present = false;
//    bool target_perc_property_present = false;
    for(const auto& obj : search_objects) {
		if(obj.is_target) {
			if(obj.perc_color == 'G') n_target_illusory_distractors++;
			else if(obj.sens_color == 'R' && obj.perc_color == '.')  n_target_illusory_blanks++;
			}
		else  { // not a target
			if(obj.perc_color == 'R') n_illusory_targets++;
			}
//        if(obj.sens_color == 'R')
//            target_sens_property_present = true;
//        if(obj.perc_color == 'R')
//            target_perc_property_present = true;
        if(target_status_trace_output) {
			if(obj.is_target && obj.sens_color == '.' && obj.perc_color == '.')
				cout << trial_g << ' ' << n_fixations << " ---CSF " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target color is neither available nor visible" << endl;
			if(obj.is_target && obj.sens_color == '.' && obj.perc_color == 'R')
				cout << trial_g << ' ' << n_fixations << " ---CSF " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target color is not available but is visible" << endl;
	// 		if(obj.is_target && obj.sens_color == 'R' && obj.perc_color != 'R')
	//            cout << trial_g << ' ' << n_fixations << " ---CSF " << obj.name << ' ' << obj.ID << 				' ' << obj.eccentricity << " target color is available but not visible" << endl;
	//       if(obj.is_target && obj.perc_color == 'R')
	//            cout << trial_g << ' ' << n_fixations << " +++CSF " << obj.name << ' ' << obj.ID << 				' ' << obj.eccentricity << " target visible" << endl;
			if(obj.is_target && obj.sens_color == 'R' && obj.perc_color == '.')
				cout << trial_g << ' ' << n_fixations << " ???CSF " << obj.name << ' ' << obj.ID <<  " target color available but blank is visible" << endl;
			if(obj.is_target && obj.sens_color == 'R' && obj.perc_color == 'G')
				cout << trial_g << ' ' << n_fixations << " ???CSF " << obj.name << ' ' << obj.ID <<  " target color available but target looks like illusory distractor" << endl;
			if(obj.is_target && obj.sens_color == '.' && obj.perc_color == 'G')
				cout << trial_g << ' ' << n_fixations << " ---CSF " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target color not available but target looks like illusory distractor" << endl;
			if(!obj.is_target && obj.perc_color == 'R')
				cout << trial_g << ' ' << n_fixations << " !!!CSF " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " illusory target visible" << endl;
			}
		}
/*	if(target_sens_property_present) {
		n_target_sens_property_present++;
		}
    if(target_perc_property_present) {
		n_target_perc_property_present++;
		}
    if(target_sens_property_present && !target_perc_property_present) {
		n_target_property_overwritten++;
		}
*/
}

void Visual_model::check_COC_target_status(const Search_objects_t& search_objects)
{
    for(const auto& obj : search_objects) {
		if(obj.is_target) {
			if ((obj.perc_color == 'G' && obj.perc_orientation == 'V') ||
				(obj.perc_color == 'R' && obj.perc_orientation == 'H')) {
					n_target_illusory_distractors++;
				}
			else if ((obj.sens_color == 'R' && obj.perc_color == '.') ||
					 (obj.sens_orientation == 'V' && obj.perc_orientation == '.')) {
					n_target_illusory_blanks++;
				}
			}
		else { // not a target
			if(obj.perc_color == 'R' && obj.perc_orientation == 'V') {
				n_illusory_targets++;
				}
			}
	   if(target_status_trace_output) {
			if(obj.is_target && obj.perc_color == 'R' && obj.perc_orientation == 'V')
				cout << trial_g << ' ' << n_fixations << " +++COC " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target visible" << endl;
			if(obj.is_target && ((obj.sens_color == 'R' && obj.perc_color != 'R') ||
								(obj.sens_orientation == 'V' && obj.perc_orientation != 'V')))
				cout << trial_g << ' ' << n_fixations << " ---COC " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target masked" << endl;
			if(!obj.is_target && obj.perc_color == 'R' && obj.perc_orientation == 'V')
				cout << trial_g << ' ' << n_fixations << " !!!COC " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " illusory target visible" << endl;
			}
		}
}

void Visual_model::check_SHP_target_status(const Search_objects_t& search_objects)
{
    bool target_masked = false;
    bool illusory_target_visible = false;
    bool all_are_distractors = true;
    for(const auto& obj : search_objects) {
		if(obj.is_target) {
			if(obj.perc_shape == '5') n_target_illusory_distractors++;
			else if(obj.sens_shape == '2' && obj.perc_shape == '.') n_target_illusory_blanks++;
			}
		else { // not a target
			if(obj.perc_shape == '2') n_illusory_targets++;
			}
		if(target_status_trace_output) {
			if(obj.perc_shape != '5')
				all_are_distractors = false;
			if(obj.is_target && obj.perc_shape == '2')
				cout << trial_g << ' ' << n_fixations << " +++SHP " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target visible" << endl;
			if(obj.is_target && obj.sens_shape == '2' && obj.perc_shape == '5') {
				target_masked = true;
				cout << trial_g << ' ' << n_fixations << " ---SHP " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target shape sensed but overwritten with distractor shape" << endl;
				}
			if(obj.is_target && obj.sens_shape == '.' && obj.perc_shape == '5') {
				target_masked = true;
				cout << trial_g << ' ' << n_fixations << " ---SHP " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " target shape unsensed but overwritten with distractor shape" << endl;
				}
			if(!obj.is_target && obj.perc_shape == '2') {
				illusory_target_visible = true;
				cout << trial_g << ' ' << n_fixations << " !!!SHP " << obj.name << ' ' << obj.ID << ' ' << obj.eccentricity << " illusory target visible" << endl;
				}
			}
		}
		if(target_status_trace_output && all_are_distractors)
        cout << trial_g << ' ' << n_fixations << " ???SHP All objects appear to be distractors" << endl;
 //   if(target_masked && !illusory_target_visible)
//        cout << trial_g << ' ' << n_fixations << " !!!SHP target feature was overwritten" << endl;

}

// Calculate the Normal CDF value for x using mean m and sd s
// using the formula in terms of the <cmath> erf function
// see Wikipedia, Normal Distribution, for this equivalence
// with the cmath>
// this gives the probability that a random deviate will be <= x
double normal_cdf(double x, double mean, double sd)
{
    return (
    0.5 * (1 + erf((x - mean)/ (sd * sqrt(2.0))))
    );
}

bool LinearConstSD_available(double eccentricity, GU::Size size, double a, double b, double sd)
{
    // use the average of horizontal and vertical size
    double obj_size = (size.h + size.v) / 2.0;
    double mean = a + b * eccentricity;
    return lapsed_gaussian_detection_function(obj_size, mean, sd, 0.); // lapse_prob = 0;
}

double LinearConstSD_detection_probability(double eccentricity, GU::Size size, double a, double b, double sd)
{
    // use the average of horizontal and vertical size
    double obj_size = (size.h + size.v) / 2.0;
    double mean = a + b * eccentricity;
    //return lapsed_gaussian_detection_function(obj_size, mean, sd, 0.); // lapse_prob = 0;
    return normal_cdf(obj_size, mean, sd);
}



------------------------------------------------------------
// Strategy.h
//  Strategy.h
//  SimpleCrowdingModel
//
//  Created by David Kieras on 4/9/18.
//

#ifndef STRATEGY_H
#define STRATEGY_H

#include "Program_constants.h"
#include "Search_object.h"
#include "Display_creator.h"
#include "Trial_statistics.h"
#include <string>
#include <set>

class Visual_model;

class Strategy {
public:
    Strategy(Display_creator& display_creator_, Condition_e cond_, int set_size_, Polarity_e polarity_, Trial_statistics& statistics_)
        : display_creator(display_creator_), condition(cond_), set_size(set_size_), polarity(polarity_), statistics(statistics_), last_response(false)
        {}
    static std::string get_strategy_id_string();
    void do_trial();
    
    


private:
    Display_creator& display_creator;
    Condition_e condition;
    int set_size;
    Polarity_e polarity;
    Trial_statistics& statistics;

    bool response = false;  // response made on the trial
    bool error = false; // true if error made on the trial
    // maintain last response information over trials
    bool last_response = false; // absent is false, present is true;
    

    Search_objects_t search_objects; // available throughout this module
    Search_objects_t::iterator nominated_object_it;
    Search_objects_t::iterator fixated_object_it;
    int fixated_object_ID;  // ID of the Search_object that is target of eye movement

    double move_eye_to_object(Visual_model& visual_model,Search_object& object_to_fixate);


    // returns motor response time
    double make_positive_response();
    double make_negative_response();
    
    // trial data collection -- must be cleared at each trial start
    int n_repeat_fixations = 0;
    int n_target_fixations = 0;
    std::set<int> fixated_objects; // keeps track of which objects have been fixated


};

#endif /* STRATEGY_H */

------------------------------------------------------------
// Point.h
// Point is a set of (x, y) coordinates. 

#ifndef POINT_H
#define POINT_H

namespace Geometry_Utilities {

/* Point */
// A Point contains an (x, y) pair to represent coordinates
struct Point
{
	double x;
	double y;

	explicit Point (double in_x = 0., double in_y = 0.) :
		x(in_x), y(in_y)
		{}

	// compare two Points
	bool operator== (const Point& rhs) const
		{return (x == rhs.x && y == rhs.y);}
 
	bool operator!= (const Point& rhs) const
		{return !(*this == rhs);}
		
	bool operator< (const Point& rhs) const
		{return (x == rhs.x) ? (y < rhs.y) : (x < rhs.x);}
		
	bool operator<= (const Point& rhs) const
		{return (*this < rhs) || (*this == rhs);}
		
	bool operator> (const Point& rhs) const
		{return !(*this <= rhs);}

	bool operator>= (const Point& rhs) const
		{return !(*this < rhs);}
};

}	// end namepsace

#endif

------------------------------------------------------------
// Ostream_format_saver.h
#ifndef OSTREAM_FORMAT_SAVER_H
#define OSTREAM_FORMAT_SAVER_H
#include <iostream>

// An RAII-concept object for saving ostream formatting states

class Ostream_format_saver {
public:
Ostream_format_saver(std::ostream& os_) :
	os(os_), old_flags(os.flags()), old_precision(os.precision()) {}
~Ostream_format_saver() 
	{
		os.flags(old_flags);
		os.precision(old_precision);
	}
private:
	std::ostream& os;
	std::ios::fmtflags old_flags;
	std::streamsize old_precision;
};
#endif

------------------------------------------------------------
// Utility_templates.h
#ifndef UTILITY_TEMPLATES_H
#define UTILITY_TEMPLATES_H

/* See Scott Meyers, Effective STL, Item 7 for a discussion
of this general purpose Function Object class - will work for deleting any pointer */
struct Delete {
template<typename T>
void operator() (const T* ptr) const
	{
		delete ptr;
	}
};

struct Delete_second {
template<typename T>
void operator() (const T thePair) const
	{
		delete thePair.second;
	}
};

#endif

------------------------------------------------------------
// View_base.h
#ifndef VIEW_BASE_H
#define VIEW_BASE_H

#include <list>
#include <string>
#include "Geometry.h"
namespace GU = Geometry_Utilities;

/* Following the Model-View-Controller pattern, View_base.h is an abstract 
base class for views used in the EPIC software. 
Each derived class will respond to a particular kind of event notification 
provided by the Model. 
*/


class Human_processor;
class Human_subprocessor;
class Symbol;
struct Speech_word;

class View_base {
public:
	View_base() :
		human_ptr(0)
		{}
	virtual ~View_base();

	void set_attached_human(Human_processor * human_ptr_)
		{human_ptr = human_ptr_;}
		
	virtual void clear() = 0; // must override

	virtual void notify_eye_movement(GU::Point)
		{}
	virtual void notify_object_appear(const Symbol&, GU::Point, GU::Size)
		{}
//	virtual void notify_object_appear(const Symbol&, GU::Point)
//		{}
	virtual void notify_object_disappear(const Symbol&)
		{}
	virtual void notify_object_reappear(const Symbol&)
		{}
	virtual void notify_erase_object(const Symbol&)
		{}
	virtual void notify_visual_location_changed(const Symbol&, GU::Point)
		{}
	virtual void notify_visual_size_changed(const Symbol&, GU::Size)
		{}
	virtual void notify_visual_property_changed(const Symbol&, const Symbol&, const Symbol&)
		{}
	virtual void notify_auditory_stream_appear(const Symbol&, double, double, GU::Point)
		{}
	virtual void notify_auditory_stream_disappear(const Symbol&)
		{}
	virtual void notify_auditory_stream_location_changed(const Symbol&, GU::Point)
		{}
	virtual void notify_auditory_stream_pitch_changed(const Symbol&, double)
		{}
	virtual void notify_auditory_stream_loudness_changed(const Symbol&, double)
		{}
	virtual void notify_auditory_stream_size_changed(const Symbol&, GU::Size)
		{}
	virtual void notify_auditory_stream_property_changed(const Symbol&, const Symbol&, const Symbol&)
		{}
	virtual void notify_auditory_sound_start(const Symbol&, const Symbol&, long, GU::Point, const Symbol&, double)
		{}
//	virtual void notify_auditory_speech_start(const Symbol&, const Symbol&, long, const Symbol&, const Symbol&, const Symbol&, double)
//		{}
	virtual void notify_auditory_speech_start(const Speech_word&)
		{}
	virtual void notify_auditory_sound_stop(const Symbol&)
		{}
	virtual void notify_erase_sound(const Symbol&)
		{}
	virtual void notify_auditory_sound_property_changed(const Symbol&, const Symbol&, const Symbol&)
		{}
	virtual void notify_append_text(const std::string&)
		{}
	virtual void notify_time(long)
		{}
	
	// attach or detach this View to or from the specified processor
	void attach_to(Human_subprocessor * processor_ptr);
	void detach_from(Human_subprocessor * processor_ptr);
	// detach this view from all current processors
	void detach_from_all();

protected:
	Human_processor * human_ptr;	// a reference to the connected model - assumed to exist when constructed
	// list of subprocessors that this view has been attached to - used by dtor
	std::list<Human_subprocessor *> processor_ptrs;
};

#endif


------------------------------------------------------------
// Exception.h
/* 
Exception is a base class for exceptions that contains a std::string
with the std::exception interface for accessing its contents. Note that 
the value returned by what() is valid only as long as the exception object 
exists.  The std::string member is protected to allow subclasses to access
it directly for more flexibility in initialization.
*/

#ifndef MY_EXCEPTION_H
#define MY_EXCEPTION_H

#include <exception>
#include <string>

class Exception : public std::exception {
public:
	Exception() : msg("Unspecified Exception")
		{}
	Exception(const std::string& msg_) :
		msg(msg_)
		{}
	~Exception() throw()
		{}
	const char * what() const noexcept override
		{return msg.c_str();}
protected:
	std::string msg;
};

#endif

------------------------------------------------------------
// Output_tee.h
#ifndef OUTPUT_TEE_H
#define OUTPUT_TEE_H

#include "View_base.h"
#include "Assert_throw.h"

#include <iostream>
#include <sstream>
#include <list>

/*

The Output_tee class allows you to split output among as many output streams
and View_bases as you wish.  It keeps a list of pointers to ostream objects and View_bases, 
and overloads operator<< to output to each stream in the list, and to an internal
ostringstring object whose contents are written a line-at-a-time to each View_base.

This class assumes the pointed-to stream objects and View_bases exist. 
If the stream object or View_base is deallocated, it should first be taken out of the list.
No check is made for putting the same stream or View_base in the list more than once. 

Note that the ostringstream object is always present, but is only written to or manipulated if
there are View_bases in the list. 
*/


class Output_tee {
public:
	// return true if the list of streams or View_bases is non-empty;
	// if the list of streams is non-empty, at least the first must be good
	// usage:
	// if(Normal_out)
	//	Normal_out << stuff;
	operator bool() const
		{
			return (
				(!stream_ptr_list.empty() && stream_ptr_list.front()->good())
			 	|| 
			 	!view_ptr_list.empty());
		}
	
	// put a stream into the list
	void add_stream(std::ostream& os)
		{
			std::ostream * os_ptr = &os;
			stream_ptr_list.push_back(os_ptr); 
		}

	// take a stream out of the list		
	void remove_stream(std::ostream& os)
		{
			stream_ptr_list.remove(&os);
		}
	// is a stream in the list
	bool is_present(std::ostream& os)
		{
			for(stream_ptr_list_t::iterator it = stream_ptr_list.begin(); it != stream_ptr_list.end(); ++it)
				if(*it == &os)
					return true;
			return false;
		}

	// put an View_base into the list
	void add_view(View_base * view_ptr)
		{
			view_ptr_list.push_back(view_ptr); 
		}

	// take an View_base out of the list		
	void remove_view(View_base * view_ptr)
		{
			view_ptr_list.remove(view_ptr);
		}

	// is a view in the list
	bool is_present(View_base * view_ptr)
		{
			for(view_ptr_list_t::iterator it = view_ptr_list.begin(); it != view_ptr_list.end(); ++it)
				if(*it == view_ptr)
					return true;
			return false;
		}

	// check if there are stream pointers in the stream pointer list - or window pointers - if so, output being produced
	bool output_enabled()
		{
			return !stream_ptr_list.empty() || !view_ptr_list.empty();
		}
	
	// templated member function applies for any type being output
	template <class T>
	Output_tee& operator<< (const T& x)
	{
		std::ostream * os_ptr;
		for (stream_ptr_list_t::iterator it = stream_ptr_list.begin(); it != stream_ptr_list.end(); it++) {
			os_ptr = *it;
			if ((*os_ptr).good())
				*os_ptr << x;
			}

		// do nothing further if no views in the list
		if(view_ptr_list.empty())
			return *this;
		
		view_line_buffer << x;

		return *this;
	}

	// specialize for manipulators defined on ios_base
	Output_tee& operator<< (std::ios_base& (*manip) (std::ios_base&))
	{
		std::ostream * os_ptr;
		for (stream_ptr_list_t::iterator it = stream_ptr_list.begin(); it != stream_ptr_list.end(); it++) {
			os_ptr = *it;
			if ((*os_ptr).good())
				manip(*os_ptr);
			}

		// do nothing further if no views in the list
		if(view_ptr_list.empty())
			return *this;

		os_ptr = &view_line_buffer;
		manip(*os_ptr);

		return *this;	
	}

	// specialize for manipulators defined on ostream
	Output_tee& operator<< (std::ostream& (*manip) (std::ostream&))
	{
		std::ostream * os_ptr;
		for (stream_ptr_list_t::iterator it = stream_ptr_list.begin(); it != stream_ptr_list.end(); it++) {
			os_ptr = *it;
			if ((*os_ptr).good())
				manip(*os_ptr);
			}
			
		// do nothing further if no views in the list
		if(view_ptr_list.empty())
			return *this;
		
		os_ptr = &view_line_buffer;
		manip(*os_ptr);
		
		// if last character was a newline, then write the line to the views and 
		// empty the stream (simulates effect of "endl")
		const std::string& s = view_line_buffer.str();
		std::size_t n = s.size();
		if(s[n - 1] != '\n')
			return *this;

		for (view_ptr_list_t::iterator it = view_ptr_list.begin(); it != view_ptr_list.end(); it++) {
			View_base * window_ptr = *it;
			window_ptr->notify_append_text(view_line_buffer.str());
			}

		view_line_buffer.str("");
			
		return *this;	
	}

	friend class Output_tee_format_saver;

private:
	typedef std::list<std::ostream *> stream_ptr_list_t;
	stream_ptr_list_t stream_ptr_list;
	typedef std::list<View_base *> view_ptr_list_t;
	view_ptr_list_t view_ptr_list;
	std::ostringstream view_line_buffer;

};

// This class provides for saving and restoring the stream formatting flags.
// Instantiating it saves the current state of the stream formatting in the Output_tee
// and then restores it when destructed. Use in any function that alters the stream formatting.
// It is assumed that the formatting for the first stream or the internal ostreamstringg 
// applies to all streams; if no streams or views are associated with the Output_tee at 
// the time of construction, the default settings are saved. The settings of the first stream
// are saved; if no streams, then the settings of the internal stringstream are saved.
// At the time of destruction, the streams or views will be set to the saved settings.
// So for consistency and clarity, the number of streams or views associated with the Output_tee 
// should not be altered during the lifetime of an Output_tee_format_saver.


// Usage example:
//	void func(Output_tee& ot)
//	{
//		Output_tee_format_saver s(ot);
//		ot << fixed << setprecision(2) << x << endl;
//	}	// state of ot restored at exit

class Output_tee_format_saver {
public:
	Output_tee_format_saver(Output_tee& output_tee_) :
		output_tee(output_tee_), old_flags(std::ios::fmtflags()), old_precision(0)
		{
			std::ostream * os_ptr = 0;
			// if there is an output stream, save the state of the first one; 
			// otherwise save the ostringstream state
			if(!output_tee.stream_ptr_list.empty())
				os_ptr = output_tee.stream_ptr_list.front();
			else if(!output_tee.view_ptr_list.empty())
				os_ptr = &output_tee.view_line_buffer;
			else
				return;
			
			old_flags = os_ptr->flags();
			old_precision = os_ptr->precision();
		}
	
	~Output_tee_format_saver()
		{
			// apply the saved formatting to all streams
			for (std::list<std::ostream *>::iterator it = output_tee.stream_ptr_list.begin(); it != output_tee.stream_ptr_list.end(); it++) {
				std::ostream * os_ptr = *it;
				os_ptr->flags(old_flags);
				os_ptr->precision(old_precision);
				}
			// apply the saved formatting to the ostringstream if views are present
			if(!output_tee.view_ptr_list.empty()) {
				output_tee.view_line_buffer.flags(old_flags);
				output_tee.view_line_buffer.precision(old_precision);
				}
		}
private:
	Output_tee& output_tee;
	std::ios::fmtflags old_flags;
	std::streamsize old_precision;
};


	
#endif

------------------------------------------------------------
// ObsData_1279_58_458.txt
438.462	435.126	427.313	427.411	406.044	408.964	415.339	417.36	637.497	705.792	788.678	927.458	538.883	543.142	587.876	652.544	847.987	1220.77	1729.59	2003.28	719.468	888.5	1161.04	1330.43
0.0274407	0.0242518	0.0225747	0.0111656	0.0379241	0.0360347	0.0313459	0.0432501	0.0409968	0.0267394	0.0203828	0.0282956	0.0294013	0.0581891	0.0552314	0.0802539	0.0214888	0.022574	0.0194678	0.0279716	0.0333564	0.0271202	0.0874717	0.164927

------------------------------------------------------------
// ObsData_48_310_129.txt
522.959	507.668	496.399	502.335	453.545	473.552	474.791	480.871	618.607	717.748	1044.87	1423.52	577.795	595.333	690.112	763.249	727.645	1071.26	1716.71	2262.38	628.676	749.101	989.609	1195.17
0.011375	0.00582524	0.00330396	0.00183824	0.003663	0.00880064	0.00993282	0.00806612	0.00395479	0.00199601	0	0	0.00416667	0.00490296	0.003663	0.00885251	0.0100976	0.00778317	0.0119291	0.00969505	0.012328	0.0104621	0.0180769	0.0333141

------------------------------------------------------------
// ObsData_All_3610_37.txt
436.697	430.042	423.620	425.694	395.417	404.149	407.709	412.718	588.905	696.319	1020.3	1390.83	547.879	579.502	663.931	756.67	1013.71	1560.34	2478.54	3001.41	817.514	1019.79	1409.28	1723.76																							
0.0181506	0.0159622	0.013423	0.00712618	0.0226247	0.0232149	0.0216471	0.030606	0.0134109	0.0139093	0.006	0.00540906	0.0160605	0.0118495	0.0197407	0.0285821	0.00612958	0.00404278	0.00738684	0.00200401	0.00895976	0.0137485	0.056776	0.0828176

------------------------------------------------------------
// ObsData_AllSubjects.txt
436.697	430.042	423.620	425.694	395.417	404.149	407.709	412.718	565.219	631.131	786.914	954.531	514.505	534.385	591.329	650.208	814.396	1191.700	1824.880	2244.170	686.874	838.747	1119.240	1331.720
0.0181506	0.0159622	0.013423	0.00712618	0.0226247	0.0232149	0.0216471	0.030606	0.0185632	0.0154087	0.0118923	0.0114084	0.0217671	0.0261225	0.0298522	0.0492985	0.0169315	0.0127917	0.0123303	0.0134132	0.0191515	0.0183662	0.0517755	0.0924684

------------------------------------------------------------
// ObsData_356_1279_37.txt
376.835	371.513	370.176	372.31	342.496	351.459	352.815	361.092	560.329	602.544	694.199	757.331	520.131	542.905	592.301	631.43	1013.71	1560.34	2478.54	3001.41	817.514	1019.79	1409.28	1723.76
0.0102809	0.0116673	0.00796687	0.00526556	0.0148665	0.0157314	0.0165249	0.0287737	0.0084841	0.00447409	0.00404513	0.000984439	0.0114116	0.00994225	0.0124725	0.0366828	0.00612958	0.00404278	0.00738684	0.00200401	0.00895976	0.0137485	0.056776	0.0828176

------------------------------------------------------------
// ObsData_1279_458_458.txt
438.462	435.126	427.313	427.411	406.044	408.964	415.339	417.36	548.053	604.06	677.148	781.169	473.631	477.908	517.431	568.784	847.987	1220.77	1729.59	2003.28	719.468	888.5	1161.04	1330.43
0.0274407	0.0242518	0.0225747	0.0111656	0.0379241	0.0360347	0.0313459	0.0432501	0.0371544	0.0314875	0.0282475	0.0313064	0.0412811	0.0619693	0.0631367	0.0868359	0.0214888	0.022574	0.0194678	0.0279716	0.0333564	0.0271202	0.0874717	0.164927

------------------------------------------------------------