ReadParams.cpp

/** \file  Params.cpp
 *  \brief Support for parsing parameter files
 */

#include "ReadParams.h"

const double ICDF_START = 100.0;


ParamMap Params::read_params_map(const char* file)
{
	std::ifstream fin(file);  // The file to be read
	std::string buf;       // Buffer for a line of text

	ParamMap param_map;   // Resulting map of (key, value)
	ParamIter iter;       // Iterator, for checking duplicate keys
	std::string key;      // Current key being processed
	std::string value;    // Current value being built
	bool skipping = true;
	while (std::getline(fin, buf))    // Each line. (buffer includes new-line)
	{
		buf.erase(0, buf.find_first_not_of(" \t\n\r"));     // Remove lead space
		buf.erase(buf.find_last_not_of(" \t\n\r") + 1);     // ... and trailing.

		if (skipping) 
		{   // Currently, we're not doing anything interesting.
			if ((buf.length() > 1) && (buf.compare(0, 1, "[") == 0)) // We found a key
			{
				skipping = false;
				key = buf.substr(1, buf.length() - 2);
			}
		}
		else
		{    // Not skipping...
			if ((buf.length() == 0) || (buf.compare(0, 1, "*") == 0) ||
				(buf.compare(0, 1, "^") == 0) || (buf.compare(0, 1, "=") == 0) || (buf.compare(0, 1, "[") == 0)) 
			{
				value.erase(0, value.find_first_not_of(" \t\n\r"));
				value.erase(value.find_last_not_of(" \t\n\r") + 1);

				// ERR if the key already exists.
				iter = param_map.find(key);
				if (iter != param_map.end()) {
					ERR_CRITICAL_FMT("Duplicate parameter values for %s\n(1):%s\n(2):%s\n",
						key.c_str(), iter->second.c_str(), value.c_str());
				}
				param_map.insert(ParamPair(key, value));
				value.clear();
				key.clear();
				if (buf.compare(0, 1, "[") == 0) 
				{
					key = buf.substr(1, buf.length() - 2);
				}
				else
				{
					skipping = true;
				}
			}
			else 
			{
				value.append(buf);
				value.append("\n");
			}
		}
	}

	if (!key.empty()) 
	{
		value.erase(0, value.find_first_not_of(" \t\n\r"));
		value.erase(value.find_last_not_of(" \t\n\r") + 1);

		iter = param_map.find(key);
		if (iter != param_map.end()) {
			ERR_CRITICAL_FMT("Duplicate parameter values for %s\n(1):%s\n(2):%s\n",
				key.c_str(), iter->second.c_str(), value.c_str());
		}

		param_map.insert(ParamPair(key, value));
	}

	fin.close();
	return param_map;
}

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

int Params::parse_int(std::string s, std::string param)
{
	try
	{
		std::string::size_type idx;
		int res = std::stoi(s, &idx);
		if (fabs(res - std::stod(s, &idx)) > 1e-10)
		{
			ERR_CRITICAL_FMT("Error - %s appears to be double, not int - value '%s'", param.c_str(), s.c_str());
		}
		return res;
	}
	catch (const std::exception& e)
	{
		ERR_CRITICAL_FMT("Error %s parsing int %s - value '%s'\n", e.what(), param.c_str(), s.c_str());
	}
}

double Params::parse_double(std::string s, std::string param)
{
	try
	{
		std::string::size_type idx;
		double result = std::stod(s, &idx);
		return result;
	}
	catch (const std::exception& e)
	{
		ERR_CRITICAL_FMT("Error %s parsing double %s - value '%s'\n", e.what(), param.c_str(), s.c_str());
	}
}

std::string Params::clp_overwrite(std::string value, Param* P) {
	if (value.at(0) != '#')
	{
		return value;
	} 
  int clp_no = Params::parse_int(value.substr(1, std::string::npos), value);
	if ((clp_no < 0) || (clp_no > 99)) 
	{
		ERR_CRITICAL_FMT("CLP %d is out of bounds reading parameters\n", clp_no);
	}
	return std::to_string(P->clP[clp_no]);
}

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

std::string Params::lookup_param(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, Param* P, bool search_clp)
{
	ParamIter iter = params.find(param_name);
	if (iter != params.end())
	{
		return search_clp ? Params::clp_overwrite(iter->second, P) : iter->second;
	}

	iter = fallback.find(param_name);
	if (iter != fallback.end()) {
		return search_clp ? Params::clp_overwrite(iter->second, P) : iter->second;
	}
	if (base == fallback)
	{
		return "NULL";
	}

	iter = base.find(param_name);
	if (iter != base.end())
	{
		return search_clp ? Params::clp_overwrite(iter->second, P) : iter->second;
	}
	return "NULL";
}


std::string Params::lookup_param(ParamMap &fallback, ParamMap &params, std::string param_name, Param* P, bool search_clp)
{
  return Params::lookup_param(fallback, fallback, params, param_name, P, search_clp);
}

std::string Params::lookup_param_clp(ParamMap& base, ParamMap& fallback, ParamMap& params, std::string param_name, Param* P)
{
	return Params::lookup_param(base, fallback, params, param_name, P, true);
}

std::string Params::lookup_param_clp(ParamMap& fallback, ParamMap& params, std::string param_name, Param* P)
{
	return Params::lookup_param(fallback, fallback, params, param_name, P, true);
}

std::string Params::lookup_param_non_clp(ParamMap& base, ParamMap& fallback, ParamMap& params, std::string param_name, Param* P)
{
	return Params::lookup_param(base, fallback, params, param_name, P, false);
}

std::string Params::lookup_param_non_clp(ParamMap& fallback, ParamMap& params, std::string param_name, Param* P)
{
	return Params::lookup_param(fallback, fallback, params, param_name, P, false);
}


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

bool Params::param_found(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name)
{

	ParamIter iter = params.find(param_name);
	if (iter != params.end())
	{
		return true;
	}

	iter = fallback.find(param_name);
	if (iter != fallback.end())
	{
		return true;
	}

	if (base == fallback)
	{
		return false;
	}

	iter = base.find(param_name);
	return (iter != base.end());
}

bool Params::param_found(ParamMap &fallback, ParamMap &params, std::string param_name)
{
	return Params::param_found(fallback, fallback, params, param_name);
}

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

double Params::get_double(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, double default_value, bool err_on_missing, Param* P)
{
	std::string str_value = Params::lookup_param_clp(base, fallback, params, param_name, P);

	if (str_value.compare("NULL") != 0)
	{
		return Params::parse_double(str_value, param_name);
	}
	if (err_on_missing)
	{
		ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
	} 
	return default_value;
}

double Params::get_double(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, double default_value, Param* P)
{
	return Params::get_double(base, fallback, params, param_name, default_value, false, P);
}

double Params::get_double(ParamMap &fallback, ParamMap &params, std::string param_name, double default_value, Param* P)
{
	return Params::get_double(fallback, fallback, params, param_name, default_value, false, P);
}

double Params::req_double(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, Param* P)
{
	return Params::get_double(base, fallback, params, param_name, 0, true, P);
}

double Params::req_double(ParamMap &fallback, ParamMap &params, std::string param_name, Param* P)
{
	return Params::req_double(fallback, fallback, params, param_name, P);
}


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

int Params::get_int(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, int default_value, bool err_on_missing, Param* P, bool force_fail)
{
	std::string str_value = Params::lookup_param_clp(base, fallback, params, param_name, P);
	
	if (!force_fail && (str_value.compare("NULL") != 0))
	{
		return Params::parse_int(str_value, param_name);
	}
	if (err_on_missing)
	{
		ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
	}
	return default_value;
}

int Params::get_int(ParamMap &fallback, ParamMap &params, std::string param_name, int default_value, Param* P)
{
	return Params::get_int(fallback, fallback, params, param_name, default_value, false, P, false);
}

int Params::get_int(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, int default_value, Param* P)
{
	return Params::get_int(base, fallback, params, param_name, default_value, false, P, false);
}

int Params::get_int_ff(bool force_fail, ParamMap& fallback, ParamMap& params, std::string param_name, int default_value, Param* P)
{
	return Params::get_int(fallback, fallback, params, param_name, default_value, false, P, force_fail);
}

int Params::req_int(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, Param* P)
{
	return Params::get_int(base, fallback, params, param_name, 0, true, P, false);
}

int Params::req_int(ParamMap &fallback, ParamMap &params, std::string param_name, Param* P)
{
	return Params::req_int(fallback, fallback, params, param_name, P);
}

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

void Params::get_double_vec(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, double* array, int expected, double default_value, int default_size, bool err_on_missing, Param* P, bool force_fail)
{
	std::string str_value = Params::lookup_param_non_clp(base, fallback, params, param_name, P);

	if ((str_value.compare("NULL") != 0) && (!force_fail))
	{
		std::stringstream str_stream(str_value);
		std::string buffer;
		int index = 0;
		while (str_stream >> buffer)
		{
			if (!buffer.empty())
			{
			  if (index < expected) array[index] = Params::parse_double(Params::clp_overwrite(buffer, P), param_name);
			  index++;
			}
		}
		if (index != expected)
		{
			Files::xfprintf_stderr("Warning - Extra elements for %s (%d - only needed %d)\n", param_name.c_str(), index, expected);
		}
		return;
	}
	if (err_on_missing)
	{
		ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
	}
	for (int i = 0; i < default_size; i++) array[i] = default_value;
}

void Params::get_double_vec(ParamMap &fallback, ParamMap &params, std::string param_name, double* array, int expected, double default_value, int default_size, Param* P)
{
	Params::get_double_vec(fallback, fallback, params, param_name, array, expected, default_value, default_size, false, P, false);
}

void Params::req_double_vec(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, double* array, int expected, Param* P)
{
	Params::get_double_vec(base, fallback, params, param_name, array, expected, 0, 0, true, P, false);
}

void Params::req_double_vec(ParamMap &fallback, ParamMap &params, std::string param_name, double* array, int expected, Param* P)
{
	Params::get_double_vec(fallback, fallback, params, param_name, array, expected, 0, 0, true, P, false);
}

void Params::get_double_vec_ff(bool force_fail, ParamMap &fallback, ParamMap &params, std::string param_name, double* array, int expected, double default_value, Param* P)
{
	Params::get_double_vec(fallback, fallback, params, param_name, array, expected, default_value, expected, false, P, force_fail);
}

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

void Params::get_int_vec(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, int default_value, int default_size, bool err_on_missing, Param* P, bool force_fail)
{
	std::string str_value = Params::lookup_param_non_clp(base, fallback, params, param_name, P);

	if ((str_value.compare("NULL") != 0) && (!force_fail))
	{
		std::stringstream str_stream(str_value);
		std::string buffer;
		int index = 0;
		while (str_stream >> buffer)
		{
			if (!buffer.empty())
			{
				int result = Params::parse_int(Params::clp_overwrite(buffer, P), param_name);
				if (index < expected)
				{
					array[index] = result;
				}
				index++;
			}
		}
		if (index != expected)
		{
			Files::xfprintf_stderr("Warning - Extra elements for %s (%d - only needed %d)\n", param_name.c_str(), index, expected);
		}
		return;
	}
	if (err_on_missing)
	{
		ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
	}
	for (int i = 0; i < default_size; i++) array[i] = default_value;
}

void Params::get_int_vec(ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, int default_value, int default_size, bool err_on_missing, Param* P, bool force_fail)
{
	Params::get_int_vec(fallback, fallback, params, param_name, array, expected, default_value, default_size, err_on_missing, P, force_fail);
}

void Params::get_int_vec(ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, int default_value, int default_size, Param* P)
{
	Params::get_int_vec(fallback, fallback, params, param_name, array, expected, default_value, default_size, false, P, false);
}

void Params::req_int_vec(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, Param* P)
{
	Params::get_int_vec(base, fallback, params, param_name, array, expected, 0, 0, true, P, false);
}

void Params::req_int_vec(ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, Param* P)
{
	Params::get_int_vec(fallback, fallback, params, param_name, array, expected, 0, 0, true, P, false);
}

void Params::get_int_vec_ff(bool force_fail, ParamMap &fallback, ParamMap &params, std::string param_name, int* array, int expected, int default_value, Param* P)
{
	Params::get_int_vec(fallback, fallback, params, param_name, array, expected, default_value, expected, false, P, force_fail);
}

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

int Params::req_string_vec(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, char** array, int expected, Param* P)
{

	std::string str_value = Params::lookup_param_non_clp(fallback, params, param_name, P);
	if (str_value.compare("NULL") != 0)
	{
		std::stringstream str_stream(str_value);
		std::string buffer;
		int index = 0;
		while (str_stream >> buffer)
		{
			if (!buffer.empty())
			{
				array[index] = new char[buffer.length() + 1];
				strcpy(array[index], Params::clp_overwrite(buffer, P).c_str());
				index++;
			}
		}
		return index;
	}
	ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
}

int Params::req_string_vec(ParamMap &fallback, ParamMap &params, std::string param_name, char** array, int expected, Param* P)
{
	return Params::req_string_vec(fallback, fallback, params, param_name, array, expected, P);
}

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

void Params::get_double_matrix(ParamMap &base, ParamMap &fallback, ParamMap &params, std::string param_name, double** array, int sizex, int sizey, double default_value, bool err_on_missing, Param* P)
{
	std::string str_value = Params::lookup_param_non_clp(base, fallback, params, param_name, P);

	if (str_value.compare("NULL") != 0)
	{
		std::stringstream str_stream(str_value);
		std::string buffer;
		int x = 0;
		int y = 0;
		int count_values = 0;
		while (str_stream >> buffer)
		{
			if (!buffer.empty())
			{
				count_values++;
				if ((y < sizey) && (x < sizex)) array[x][y] = Params::parse_double(Params::clp_overwrite(buffer, P), param_name);
				x++;
				if (x == sizex)
				{
					y++;
					x = 0;
				}
			}
		}
		if (count_values != (sizex * sizey))
		{
			Files::xfprintf_stderr("Warning: Expected %d values for matrix %s - actually available; %d\n", sizex * sizey, param_name.c_str(), count_values);
		}
		return;
	}

	if (err_on_missing)
	{
		ERR_CRITICAL_FMT("Required Parameter %s not found\n", param_name.c_str());
		return;
	}

	for (int x = 0; x < sizex; x++)
	{
		for (int y = 0; y < sizey; y++)
		{
			array[x][y] = default_value;
		}
	}
}

void Params::get_double_matrix(ParamMap &fallback, ParamMap &params, std::string param_name, double** array, int sizex, int sizey, double default_value, Param* P)
{
	Params::get_double_matrix(fallback, fallback, params, param_name, array, sizex, sizey, default_value, false, P);
}

void Params::get_inverse_cdf(ParamMap fallback, ParamMap params, const char* icdf_name, InverseCdf* inverseCdf, Param* P, double start_value)
{
	Params::get_double_vec(fallback, params, icdf_name, inverseCdf->get_values(), CDF_RES + 1, 0, CDF_RES + 1, P);
	if (!Params::param_found(fallback, params, icdf_name))
	{
		inverseCdf->set_neg_log(start_value);
	}
	inverseCdf->assign_exponent();
}

double** create_2d_double(int sizex, int sizey)
{
	double** arr = new double* [sizex]();
	for (int i = 0; i < sizex; i++)	arr[i] = new double[sizey]();
	return arr;
}

double*** create_3d_double(int sizex, int sizey, int sizez)
{
	double*** arr = new double** [sizex]();
	for (int i = 0; i < sizex; i++)	arr[i] = create_2d_double(sizey, sizez);
	return arr;
}

void Params::alloc_params(Param* P)
{
	P->LocationInitialInfection				= create_2d_double(MAX_NUM_SEED_LOCATIONS, 2);
	P->WAIFW_Matrix							= create_2d_double(NUM_AGE_GROUPS, NUM_AGE_GROUPS);
	P->WAIFW_Matrix_SpatialOnly				= create_2d_double(NUM_AGE_GROUPS, NUM_AGE_GROUPS);
	P->SD_PlaceEffects_OverTime				= create_2d_double(MAX_NUM_INTERVENTION_CHANGE_TIMES, MAX_NUM_PLACE_TYPES);
	P->Enhanced_SD_PlaceEffects_OverTime	= create_2d_double(MAX_NUM_INTERVENTION_CHANGE_TIMES, MAX_NUM_PLACE_TYPES);
	P->HQ_PlaceEffects_OverTime				= create_2d_double(MAX_NUM_INTERVENTION_CHANGE_TIMES, MAX_NUM_PLACE_TYPES);
	P->PC_PlaceEffects_OverTime				= create_2d_double(MAX_NUM_INTERVENTION_CHANGE_TIMES, MAX_NUM_PLACE_TYPES);
	P->PC_PropAttending_OverTime			= create_2d_double(MAX_NUM_INTERVENTION_CHANGE_TIMES, MAX_NUM_PLACE_TYPES);
	P->HouseholdSizeDistrib					= create_2d_double(MAX_ADUNITS, MAX_HOUSEHOLD_SIZE);
	P->PropAgeGroup							= create_2d_double(MAX_ADUNITS, NUM_AGE_GROUPS);
	P->PopByAdunit							= create_2d_double(MAX_ADUNITS, 2);
	P->InvLifeExpecDist						= create_2d_double(MAX_ADUNITS, 1001);
}
/**************************************************************************************************************/

void Params::output_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->OutputAge = Params::get_int(params, pre_params, "OutputAge", 1, P);
	P->OutputSeverity = Params::get_int(params, pre_params, "OutputSeverity", 1, P);
	P->OutputSeverityAdminUnit = Params::get_int(params, pre_params, "OutputSeverityAdminUnit", 1, P);
	P->OutputSeverityAge = Params::get_int(params, pre_params, "OutputSeverityAge", 1, P);
	P->OutputAdUnitAge = Params::get_int(params, pre_params, "OutputAdUnitAge", 0, P);
	P->OutputR0 = Params::get_int(params, pre_params, "OutputR0", 0, P);
	P->OutputControls = Params::get_int(params, pre_params, "OutputControls", 0, P);
	P->OutputCountry = Params::get_int(params, pre_params, "OutputCountry", 0, P);
	P->OutputAdUnitVar = Params::get_int(params, pre_params, "OutputAdUnitVar", 0, P);
	P->OutputHousehold = Params::get_int(params, pre_params, "OutputHousehold", 0, P);
	P->OutputInfType = Params::get_int(params, pre_params, "OutputInfType", 0, P);
	P->OutputNonSeverity = Params::get_int(params, pre_params, "OutputNonSeverity", 0, P);
	P->OutputNonSummaryResults = Params::get_int(params, pre_params, "OutputNonSummaryResults", 0, P);
}

void Params::household_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	if (P->DoHouseholds == 0)
	{
		P->HouseholdTrans = 0.0;
		P->HouseholdTransPow = 1.0;
		P->HouseholdSizeDistrib[0][0] = 1.0;
		for (int i = 1; i < MAX_HOUSEHOLD_SIZE; i++)
		P->HouseholdSizeDistrib[0][i] = 0;
		return;
	}

	Params::req_double_vec(pre_params, adm_params, "Household size distribution", P->HouseholdSizeDistrib[0], MAX_HOUSEHOLD_SIZE, P);
	P->HouseholdTrans = Params::req_double(params, pre_params, "Household attack rate", P);
	P->HouseholdTransPow = Params::req_double(params, pre_params, "Household transmission denominator power", P);
	P->DoCorrectAgeDist = Params::get_int(pre_params, adm_params, "Correct age distribution after household allocation to exactly match specified demography", 0, P);
	if (P->FitIter != 0)
	{
		return;
	}
	for (int i = 1; i < MAX_HOUSEHOLD_SIZE; i++)
		P->HouseholdSizeDistrib[0][i] = P->HouseholdSizeDistrib[0][i] + P->HouseholdSizeDistrib[0][i - 1];
	P->HouseholdDenomLookup[0] = 1.0;
	for (int i = 1; i < MAX_HOUSEHOLD_SIZE; i++)
		P->HouseholdDenomLookup[i] = 1 / pow(((double)(INT64_C(1) + i)), P->HouseholdTransPow);
}

void Params::waifw_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	//if (!GetInputParameter2(params, pre_params, "WAIFW matrix", "%lf", (void*)P->WAIFW_Matrix, NUM_AGE_GROUPS, NUM_AGE_GROUPS, 0))
	if (!Params::param_found(params, pre_params, "WAIFW matrix"))
	{
		for (int i = 0; i < NUM_AGE_GROUPS; i++)
			for (int j = 0; j < NUM_AGE_GROUPS; j++)
				P->WAIFW_Matrix[i][j] = 1.0;
	}
	else
	{
		Params::get_double_matrix(params, pre_params, "WAIFW matrix", P->WAIFW_Matrix, NUM_AGE_GROUPS, NUM_AGE_GROUPS, 1.0, P);

		/* WAIFW matrix needs to be scaled to have max value of 1.
		1st index of matrix specifies host being infected, second the infector.
		Overall age variation in infectiousness/contact rates/susceptibility should be factored
		out of WAIFW_matrix and put in Age dep infectiousness/susceptibility for efficiency. */
		double t = 0;
		for (int i = 0; i < NUM_AGE_GROUPS; i++)
			for (int j = 0; j < NUM_AGE_GROUPS; j++)
				if (P->WAIFW_Matrix[i][j] > t) t = P->WAIFW_Matrix[i][j];
		if (t > 0)
		{
			for (int i = 0; i < NUM_AGE_GROUPS; i++)
				for (int j = 0; j < NUM_AGE_GROUPS; j++)
					P->WAIFW_Matrix[i][j] /= t;
		}
		else
		{
			for (int i = 0; i < NUM_AGE_GROUPS; i++)
				for (int j = 0; j < NUM_AGE_GROUPS; j++)
					P->WAIFW_Matrix[i][j] = 1.0;
		}
	}
	if (!Params::param_found(params, pre_params, "WAIFW matrix spatial infections only"))
	{
		for (int i = 0; i < NUM_AGE_GROUPS; i++)
			for (int j = 0; j < NUM_AGE_GROUPS; j++)
				P->WAIFW_Matrix_SpatialOnly[i][j] = 1.0;

		P->Got_WAIFW_Matrix_Spatial = 0;
	}
	else
	{
		Params::get_double_matrix(params, pre_params, "WAIFW matrix spatial infections only", P->WAIFW_Matrix_SpatialOnly, NUM_AGE_GROUPS, NUM_AGE_GROUPS, 0, P);
		P->Got_WAIFW_Matrix_Spatial = 1;
		/* WAIFW matrix needs to be scaled to have max value of 1.
		1st index of matrix specifies host being infected, second the infector.
		Overall age variation in infectiousness/contact rates/susceptibility should be factored
		out of WAIFW_matrix and put in Age dep infectiousness/susceptibility for efficiency. */

		double Maximum = 0;
		for (int i = 0; i < NUM_AGE_GROUPS; i++)
			for (int j = 0; j < NUM_AGE_GROUPS; j++)
				if (P->WAIFW_Matrix_SpatialOnly[i][j] > Maximum) Maximum = P->WAIFW_Matrix_SpatialOnly[i][j];
		if (Maximum > 0)
		{
			for (int i = 0; i < NUM_AGE_GROUPS; i++)
				for (int j = 0; j < NUM_AGE_GROUPS; j++)
					P->WAIFW_Matrix_SpatialOnly[i][j] /= Maximum;
		}
		else
		{
			for (int i = 0; i < NUM_AGE_GROUPS; i++)
				for (int j = 0; j < NUM_AGE_GROUPS; j++)
					P->WAIFW_Matrix_SpatialOnly[i][j] = 1.0;
		}
	}
}


///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
///// **** AIRPORT PARAMETERS
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****

void Params::airport_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->DoAirports = Params::get_int(params, pre_params, "Include air travel", 0, P);
	if (P->DoAirports == 0)  // Airports disabled => all places are not to do with airports, and we have no hotels
	{
		P->PlaceTypeNoAirNum = P->NumPlaceTypes;
		P->HotelPlaceType = P->NumPlaceTypes;
		return;
	}
	
	// When airports are activated we must have at least one airport place
	// // and a hotel type.
	P->PlaceTypeNoAirNum = Params::req_int(pre_params, adm_params, "Number of non-airport places", P);
	P->HotelPlaceType = Params::req_int(pre_params, adm_params, "Hotel place type", P);
	if (P->PlaceTypeNoAirNum >= P->NumPlaceTypes)
	{
		ERR_CRITICAL_FMT("[Number of non-airport places] parameter (%d) is greater than number of places (%d).\n", P->PlaceTypeNoAirNum, P->NumPlaceTypes);
	}
	if (P->HotelPlaceType < P->PlaceTypeNoAirNum || P->HotelPlaceType >= P->NumPlaceTypes)
	{
		ERR_CRITICAL_FMT("[Hotel place type] parameter (%d) not in the range [%d, %d)\n", P->HotelPlaceType, P->PlaceTypeNoAirNum, P->NumPlaceTypes);
	}

	P->AirportTrafficScale = Params::get_double(params, pre_params, "Scaling factor for input file to convert to daily traffic", 1.0, P);
	P->HotelPropLocal = Params::get_double(params, pre_params, "Proportion of hotel attendees who are local", 0, P);

	// If params are not specified, get_double_vec here fills with zeroes.

	Params::get_double_vec(params, pre_params, "Distribution of duration of air journeys", P->JourneyDurationDistrib, MAX_TRAVEL_TIME, 0, MAX_TRAVEL_TIME, P);
	if (!Params::param_found(params, pre_params, "Distribution of duration of air journeys"))
	{
		P->JourneyDurationDistrib[0] = 1;
	}
	Params::get_double_vec(params, pre_params, "Distribution of duration of local journeys", P->LocalJourneyDurationDistrib, MAX_TRAVEL_TIME, 0, MAX_TRAVEL_TIME, P);
	if (!Params::param_found(params, pre_params, "Distribution of duration of local journeys"))
	{
		P->LocalJourneyDurationDistrib[0] = 1;
	}
	P->MeanJourneyTime = 0;
	P->MeanLocalJourneyTime = 0;
	for (int i = 0; i < MAX_TRAVEL_TIME; i++)
	{
		P->MeanJourneyTime += ((double)(i)) * P->JourneyDurationDistrib[i];
		P->MeanLocalJourneyTime += ((double)(i)) * P->LocalJourneyDurationDistrib[i];
	}
	Files::xfprintf_stderr("Mean duration of local journeys = %lf days\n", P->MeanLocalJourneyTime);
	for (int i = 1; i < MAX_TRAVEL_TIME; i++)
	{
		P->JourneyDurationDistrib[i] += P->JourneyDurationDistrib[i - 1];
		P->LocalJourneyDurationDistrib[i] += P->LocalJourneyDurationDistrib[i - 1];
	}
	int j1 = 0;
	int j2 = 0;
	for (int i = 0; i <= 1024; i++)
	{
		double s = ((double) i) / 1024;
		while (P->JourneyDurationDistrib[j1] < s) j1++;
		P->InvJourneyDurationDistrib[i] = j1;
		while (P->LocalJourneyDurationDistrib[j2] < s) j2++;
		P->InvLocalJourneyDurationDistrib[i] = j2;
	}
}

void Params::serology_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->SeroConvMaxSens = Params::get_double(params, pre_params, "Maximum sensitivity of serology assay", 1.0, P);
	P->SeroConvP1 = Params::get_double(params, pre_params, "Seroconversion model parameter 1", 14.0, P);
	P->SeroConvP2 = Params::get_double(params, pre_params, "Seroconversion model parameter 2", 3.0, P);
	P->SeroConvSpec = Params::get_double(params, pre_params, "Specificity of serology assay", 1.0, P);
	P->InfPrevSurveyScale = Params::get_double(params, pre_params, "Scaling of modelled infection prevalence to match surveys", 1.0, P);
}


///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
///// **** SEVERITY PARAMETERS
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****

void Params::severity_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->DoSeverity = Params::get_int(params, pre_params, "Do Severity Analysis", 0, P);
	if (P->DoSeverity == 0)
	{
	  return;
	}
	P->ScaleSymptProportions = Params::get_double(params, pre_params, "Factor to scale IFR", 1.0, P);
	//// Means for icdf's.
	P->Mean_TimeToTest = Params::get_double(params, pre_params, "MeanTimeToTest", 0.0, P);
	P->Mean_TimeToTestOffset = Params::get_double(params, pre_params, "MeanTimeToTestOffset", 1.0, P);
	P->Mean_TimeToTestCriticalOffset = Params::get_double(params, pre_params, "MeanTimeToTestCriticalOffset", 1.0, P);
	P->Mean_TimeToTestCritRecovOffset = Params::get_double(params, pre_params, "MeanTimeToTestCritRecovOffset", 1.0, P);
	if (Params::get_int(params, pre_params, "Age dependent severity delays", 0, P) == 0)
	{
		P->Mean_MildToRecovery[0] = Params::req_double(params, pre_params, "Mean_MildToRecovery", P);
		P->Mean_ILIToRecovery[0] = Params::req_double(params, pre_params, "Mean_ILIToRecovery", P);
		P->Mean_SARIToRecovery[0] = Params::req_double(params, pre_params, "Mean_SARIToRecovery", P);
		P->Mean_CriticalToCritRecov[0] = Params::req_double(params, pre_params, "Mean_CriticalToCritRecov", P);
		P->Mean_CritRecovToRecov[0] = Params::req_double(params, pre_params, "Mean_CritRecovToRecov", P);
		P->Mean_ILIToSARI[0] = Params::req_double(params, pre_params, "Mean_ILIToSARI", P);
		P->Mean_ILIToDeath[0] = Params::get_double(params, pre_params, "Mean_ILIToDeath", 7.0, P);
		P->Mean_SARIToCritical[0] = Params::req_double(params, pre_params, "Mean_SARIToCritical", P);
		P->Mean_SARIToDeath[0] = Params::req_double(params, pre_params, "Mean_SARIToDeath", P);
		P->Mean_CriticalToDeath[0] = Params::req_double(params, pre_params, "Mean_CriticalToDeath", P);
		for (int AgeGroup = 1; AgeGroup < NUM_AGE_GROUPS; AgeGroup++)
		{
			P->Mean_MildToRecovery[AgeGroup] = P->Mean_MildToRecovery[0];
			P->Mean_ILIToRecovery[AgeGroup] = P->Mean_ILIToRecovery[0];
			P->Mean_SARIToRecovery[AgeGroup] = P->Mean_SARIToRecovery[0];
			P->Mean_CriticalToCritRecov[AgeGroup] = P->Mean_CriticalToCritRecov[0];
			P->Mean_CritRecovToRecov[AgeGroup] = P->Mean_CritRecovToRecov[0];
			P->Mean_ILIToSARI[AgeGroup] = P->Mean_ILIToSARI[0];
			P->Mean_ILIToDeath[AgeGroup] = P->Mean_ILIToDeath[0];
			P->Mean_SARIToCritical[AgeGroup] = P->Mean_SARIToCritical[0];
			P->Mean_SARIToDeath[AgeGroup] = P->Mean_SARIToDeath[0];
			P->Mean_CriticalToDeath[AgeGroup] = P->Mean_CriticalToDeath[0];
		}
	}
	else
	{
		Params::req_double_vec(params, pre_params, "Mean_MildToRecovery", P->Mean_MildToRecovery, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_ILIToRecovery", P->Mean_ILIToRecovery, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_SARIToRecovery", P->Mean_SARIToRecovery, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_CriticalToCritRecov", P->Mean_CriticalToCritRecov, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_CritRecovToRecov", P->Mean_CritRecovToRecov, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_ILIToSARI", P->Mean_ILIToSARI, NUM_AGE_GROUPS, P);
		Params::get_double_vec(params, pre_params, "Mean_ILIToDeath", P->Mean_ILIToDeath, NUM_AGE_GROUPS, 7.0, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_SARIToCritical", P->Mean_SARIToCritical, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_SARIToDeath", P->Mean_SARIToDeath, NUM_AGE_GROUPS, P);
		Params::req_double_vec(params, pre_params, "Mean_CriticalToDeath", P->Mean_CriticalToDeath, NUM_AGE_GROUPS, P);
	}

	//// Get InverseCDFs
	Params::get_inverse_cdf(params, pre_params, "MildToRecovery_icdf", &P->MildToRecovery_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "ILIToRecovery_icdf", &P->ILIToRecovery_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "ILIToDeath_icdf", &P->ILIToDeath_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "SARIToRecovery_icdf", &P->SARIToRecovery_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "CriticalToCritRecov_icdf", &P->CriticalToCritRecov_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "CritRecovToRecov_icdf", &P->CritRecovToRecov_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "ILIToSARI_icdf", &P->ILIToSARI_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "SARIToCritical_icdf", &P->SARIToCritical_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "SARIToDeath_icdf", &P->SARIToDeath_icdf, P, ICDF_START);
	Params::get_inverse_cdf(params, pre_params, "CriticalToDeath_icdf", &P->CriticalToDeath_icdf, P, ICDF_START);

	// If you decide to decompose Critical -> Death transition into Critical -> Stepdown and Stepdown -> Death, use the block below.
	P->IncludeStepDownToDeath = Params::get_int(params, pre_params, "IncludeStepDownToDeath", 0, P);
	if (P->IncludeStepDownToDeath == 0) /// for backwards compatibility. If Stepdown to death not included (or if unspecified), set stepdown->death = stepdown->recovery.
	{
		for (int quantile = 0; quantile <= CDF_RES; quantile++)
			P->StepdownToDeath_icdf[quantile] = P->CritRecovToRecov_icdf[quantile];
		for (int AgeGroup = 0; AgeGroup < NUM_AGE_GROUPS; AgeGroup++)
			P->Mean_StepdownToDeath[AgeGroup] = P->Mean_CritRecovToRecov[AgeGroup];
	}
	else
	{
		Params::req_double_vec(params, pre_params, "Mean_StepdownToDeath", P->Mean_StepdownToDeath, NUM_AGE_GROUPS, P);
		Params::get_inverse_cdf(params, pre_params, "StepdownToDeath_icdf", &P->StepdownToDeath_icdf, P, ICDF_START);
	}

	Params::get_double_vec(params, pre_params, "Prop_Mild_ByAge", P->Prop_Mild_ByAge, NUM_AGE_GROUPS, 0.5, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "Prop_ILI_ByAge", P->Prop_ILI_ByAge, NUM_AGE_GROUPS, 0.3, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "Prop_SARI_ByAge", P->Prop_SARI_ByAge, NUM_AGE_GROUPS, 0.15, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "Prop_Critical_ByAge", P->Prop_Critical_ByAge, NUM_AGE_GROUPS, 0.05, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "CFR_SARI_ByAge", P->CFR_SARI_ByAge, NUM_AGE_GROUPS, 0.5, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "CFR_Critical_ByAge", P->CFR_Critical_ByAge, NUM_AGE_GROUPS, 0.5, NUM_AGE_GROUPS, P);
	Params::get_double_vec(params, pre_params, "CFR_ILI_ByAge", P->CFR_ILI_ByAge, NUM_AGE_GROUPS, 0, NUM_AGE_GROUPS, P);

	//Add param to allow severity to be uniformly scaled up or down.
	for (int i = 0; i < NUM_AGE_GROUPS; i++)
	{
		P->Prop_SARI_ByAge[i] *= P->ScaleSymptProportions;
		P->Prop_Critical_ByAge[i] *= P->ScaleSymptProportions;
		P->Prop_ILI_ByAge[i] = 1.0 - P->Prop_Mild_ByAge[i] - P->Prop_SARI_ByAge[i] - P->Prop_Critical_ByAge[i];
	}
}


///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
///// **** VACCINATION PARAMETERS
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****

void Params::vaccination_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->VaccCellIncThresh = Params::get_double(params, pre_params, "Vaccination trigger incidence per cell", 1000000000, P);
	P->VaccSuscDrop = Params::get_double(params, pre_params, "Relative susceptibility of vaccinated individual", 1, P);
	P->VaccSuscDrop2 = Params::get_double(params, pre_params, "Relative susceptibility of individual vaccinated after switch time", 1, P);
	P->VaccTimeEfficacySwitch = Params::get_double(params, pre_params, "Switch time at which vaccine efficacy increases", USHRT_MAX / P->TimeStepsPerDay, P);
	P->VaccEfficacyDecay = Params::get_double(params, pre_params, "Decay rate of vaccine efficacy (per year)", 0, P);
	P->VaccEfficacyDecay /= DAYS_PER_YEAR;
	P->VaccInfDrop = Params::get_double(params, pre_params, "Relative infectiousness of vaccinated individual", 1, P);
	P->VaccMortDrop = Params::get_double(params, pre_params, "Proportion of symptomatic cases resulting in death prevented by vaccination", 0, P);
	P->VaccSympDrop = Params::get_double(params, pre_params, "Proportion of symptomatic cases prevented by vaccination", 0, P);
	P->VaccDelayMean = Params::get_double(params, pre_params, "Delay to vaccinate", 0, P);

	P->VaccTimeToEfficacy = Params::get_double(params, pre_params, "Delay from vaccination to full protection", 0, P);

	P->VaccCampaignInterval = Params::get_double(params, pre_params, "Years between rounds of vaccination", 1e10, P);
	P->VaccDosePerDay = Params::get_int(params, pre_params, "Max vaccine doses per day", -1, P);
	P->VaccCampaignInterval *= DAYS_PER_YEAR;
	P->VaccMaxRounds = Params::get_int(params, pre_params, "Maximum number of rounds of vaccination", 1, P);
	if (P->DoHouseholds != 0)
	{
		P->VaccPropCaseHouseholds = Params::get_double(params, pre_params, "Proportion of households of cases vaccinated", 0, P);
		P->VaccHouseholdsDuration = Params::get_double(params, pre_params, "Duration of household vaccination policy", USHRT_MAX / P->TimeStepsPerDay, P);
	}

	P->VaccTimeStartBase = Params::get_double(params, pre_params, "Vaccination start time", USHRT_MAX / P->TimeStepsPerDay, P);
	P->VaccProp = Params::get_double(params, pre_params, "Proportion of population vaccinated", 0, P);
	P->VaccCoverageIncreasePeriod = Params::get_double(params, pre_params, "Time taken to reach max vaccination coverage (in years)", 0, P);
	P->VaccCoverageIncreasePeriod *= DAYS_PER_YEAR;
	P->VaccTimeStartGeo = Params::get_double(params, pre_params, "Time to start geographic vaccination", 1e10, P);
	P->VaccRadius = Params::get_double(params, pre_params, "Vaccination radius", 0, P);
	P->VaccMinRadius = Params::get_double(params, pre_params, "Minimum radius from case to vaccinate", 0, P);
	P->VaccMaxCoursesBase = Params::get_double(params, pre_params, "Maximum number of vaccine courses available", 1e20, P);
	P->VaccNewCoursesStartTime = Params::get_double(params, pre_params, "Start time of additional vaccine production", USHRT_MAX / P->TimeStepsPerDay, P);
	P->VaccNewCoursesEndTime = Params::get_double(params, pre_params, "End time of additional vaccine production", USHRT_MAX / P->TimeStepsPerDay, P);
	P->VaccNewCoursesRate = Params::get_double(params, pre_params, "Rate of additional vaccine production (courses per day)", 0, P);
	P->DoMassVacc = Params::get_int(params, pre_params, "Apply mass rather than reactive vaccination", 0, P);
	if (Params::param_found(params, pre_params, "Priority age range for mass vaccination")) {
		Params::req_int_vec(params, pre_params, "Priority age range for mass vaccination", P->VaccPriorityGroupAge, 2, P);
	}
	else {
		P->VaccPriorityGroupAge[0] = 1; P->VaccPriorityGroupAge[1] = 0;
	}

	if (P->DoAdUnits == 0)
	{
		P->VaccAdminUnitDivisor = 1;
		P->VaccByAdminUnit = 0;
		return;
	}
	P->VaccByAdminUnit = Params::get_int(params, pre_params, "Vaccinate administrative units rather than rings", 0, P);
	P->VaccAdminUnitDivisor = Params::get_int(params, pre_params, "Administrative unit divisor for vaccination", 1, P);
	if ((P->VaccAdminUnitDivisor == 0) || (P->VaccByAdminUnit == 0)) P->VaccAdminUnitDivisor = 1;
}

///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
///// **** TREATMENT PARAMETERS
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****


void Params::treatment_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->DoPlaceGroupTreat = Params::get_int(params, pre_params, "Only treat mixing groups within places", 0, P);

	P->TreatCellIncThresh = Params::get_double(params, pre_params, "Treatment trigger incidence per cell", INT32_MAX, P);
	P->CaseIsolation_CellIncThresh = Params::get_double(params, pre_params, "Case isolation trigger incidence per cell", P->TreatCellIncThresh, P);
	P->HHQuar_CellIncThresh = Params::get_double(params, pre_params, "Household quarantine trigger incidence per cell", P->TreatCellIncThresh, P);

	P->TreatSuscDrop = Params::get_double(params, pre_params, "Relative susceptibility of treated individual", 1, P);
	P->TreatInfDrop = Params::get_double(params, pre_params, "Relative infectiousness of treated individual", 1, P);
	P->TreatDeathDrop = Params::get_double(params, pre_params, "Proportion of symptomatic cases resulting in death prevented by treatment", 0, P);
	P->TreatSympDrop = Params::get_double(params, pre_params, "Proportion of symptomatic cases prevented by treatment", 0, P);
	P->TreatDelayMean = Params::get_double(params, pre_params, "Delay to treat cell", 0, P);
	P->TreatCaseCourseLength = Params::get_double(params, pre_params, "Duration of course of treatment", 5, P);
	P->TreatProphCourseLength = Params::get_double(params, pre_params, "Duration of course of prophylaxis", 10, P);
	P->TreatPropCases = Params::get_double(params, pre_params, "Proportion of detected cases treated", 1, P);
	if (P->DoHouseholds != 0)
	{
		P->TreatPropCaseHouseholds = Params::get_double(params, pre_params, "Proportion of households of cases treated", 0, P);
		P->TreatHouseholdsDuration = Params::get_double(params, pre_params, "Duration of household prophylaxis policy", USHRT_MAX / P->TimeStepsPerDay, P);
	}
	// Check below - "Proportional treated" will always be ignored.
	//if (!GetInputParameter2(params, pre_params, "Proportion treated", "%lf", (void*) & (P->TreatPropRadial), 1, 1, 0)) P->TreatPropRadial = 1.0;
	//if (!GetInputParameter2(params, pre_params, "Proportion treated in radial prophylaxis", "%lf", (void*) & (P->TreatPropRadial), 1, 1, 0)) P->TreatPropRadial = 1.0;

	P->TreatPropRadial = Params::get_double(params, pre_params, "Proportion treated in radial prophylaxis", 1.0, P);
	P->TreatRadius = Params::get_double(params, pre_params, "Treatment radius", 0, P);
	P->TreatPlaceGeogDuration = Params::get_double(params, pre_params, "Duration of place/geographic prophylaxis policy", USHRT_MAX / P->TimeStepsPerDay, P);
	P->TreatTimeStartBase = Params::get_double(params, pre_params, "Treatment start time", USHRT_MAX / P->TimeStepsPerDay, P);
	if (P->DoPlaces != 0)
	{
		Params::get_double_vec(params, pre_params, "Proportion of places treated after case detected", P->TreatPlaceProbCaseId, P->NumPlaceTypes, 0, MAX_NUM_PLACE_TYPES, P);
		Params::get_double_vec(params, pre_params, "Proportion of people treated in targeted places", P->TreatPlaceTotalProp, P->NumPlaceTypes, 0, MAX_NUM_PLACE_TYPES, P);
	}
	P->TreatMaxCoursesBase = Params::get_double(params, pre_params, "Maximum number of doses available", 1e20, P);
	P->TreatNewCoursesStartTime = Params::get_double(params, pre_params, "Start time of additional treatment production", USHRT_MAX / P->TimeStepsPerDay, P);
	P->TreatNewCoursesRate = Params::get_double(params, pre_params, "Rate of additional treatment production (courses per day)", 0, P);
	P->TreatMaxCoursesPerCase = Params::get_int(params, pre_params, "Maximum number of people targeted with radial prophylaxis per case", INT32_MAX, P);

	if (P->DoAdUnits == 0)
	{
		P->TreatAdminUnitDivisor = 1;
		P->TreatByAdminUnit = 0; return;
	}
	P->TreatByAdminUnit = Params::get_int(params, pre_params, "Treat administrative units rather than rings", 0, P);
	P->TreatAdminUnitDivisor = Params::get_int(params, pre_params, "Administrative unit divisor for treatment", 1, P);
	if ((P->TreatAdminUnitDivisor == 0) || (P->TreatByAdminUnit == 0))
	{
		P->TreatByAdminUnit = 0;
		P->TreatAdminUnitDivisor = 1;
	}
}

void Params::carehome_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	P->CareHomeResidentHouseholdScaling = Params::get_double(pre_params, adm_params, "Scaling of household contacts for care home residents", 1.0, P);
	P->CareHomeResidentSpatialScaling = Params::get_double(pre_params, adm_params, "Scaling of spatial contacts for care home residents", 1.0, P);
	P->CareHomeResidentPlaceScaling = Params::get_double(pre_params, adm_params, "Scaling of between group (home) contacts for care home residents", 1.0, P);
	P->CareHomeWorkerGroupScaling = Params::get_double(pre_params, adm_params, "Scaling of within group (home) contacts for care home workers", 1.0, P);
	P->CareHomeRelProbHosp = Params::get_double(pre_params, adm_params, "Relative probability that care home residents are hospitalised", 1.0, P);
	if (P->FitIter !=0)
	{
	  return;
	}

	if (P->NumPlaceTypes > MAX_NUM_PLACE_TYPES) ERR_CRITICAL("Too many place types\n");
	P->CareHomePlaceType = Params::get_int(pre_params, adm_params, "Place type number for care homes", -1, P);
	P->CareHomeAllowInitialInfections = Params::get_int(pre_params, adm_params, "Allow initial infections to be in care homes", 0, P);
	P->CareHomeResidentMinimumAge = Params::get_int(pre_params, adm_params, "Minimum age of care home residents", 1000, P);
}

void Params::place_type_params(ParamMap adm_params, ParamMap pre_params, ParamMap params, Param* P)
{
	if (P->DoPlaces != 0)
	{
		Params::carehome_params(adm_params, pre_params, params, P);
		Params::req_double_vec(params, pre_params, "Proportion of between group place links", P->PlaceTypePropBetweenGroupLinks, P->NumPlaceTypes, P);
		Params::req_double_vec(params, pre_params, "Relative transmission rates for place types", P->PlaceTypeTrans, P->NumPlaceTypes, P);

		if (P->FitIter != 0)
		{
			return;
		}
		Params::req_int_vec(pre_params, adm_params, "Minimum age for age group 1 in place types", P->PlaceTypeAgeMin, P->NumPlaceTypes, P);
		Params::req_int_vec(pre_params, adm_params, "Maximum age for age group 1 in place types", P->PlaceTypeAgeMax, P->NumPlaceTypes, P);
		Params::req_double_vec(pre_params, adm_params, "Proportion of age group 1 in place types", P->PlaceTypePropAgeGroup, P->NumPlaceTypes, P);

		if (!Params::param_found(pre_params, adm_params, "Proportion of age group 2 in place types"))
		{
			for (int i = 0; i < MAX_NUM_PLACE_TYPES; i++)
			{
				P->PlaceTypePropAgeGroup2[i] = 0;
				P->PlaceTypeAgeMin2[i] = 0;
				P->PlaceTypeAgeMax2[i] = 1000;
			}
		}
		else
		{
			Params::req_double_vec(pre_params, adm_params, "Proportion of age group 2 in place types", P->PlaceTypePropAgeGroup2, P->NumPlaceTypes, P);
			Params::req_int_vec(pre_params, adm_params, "Minimum age for age group 2 in place types", P->PlaceTypeAgeMin2, P->NumPlaceTypes, P);
			Params::req_int_vec(pre_params, adm_params, "Maximum age for age group 2 in place types", P->PlaceTypeAgeMax2, P->NumPlaceTypes, P);
		}
		if (!Params::param_found(pre_params, adm_params, "Proportion of age group 3 in place types"))
		{
			for (int i = 0; i < MAX_NUM_PLACE_TYPES; i++)
			{
				P->PlaceTypePropAgeGroup3[i] = 0;
				P->PlaceTypeAgeMin3[i] = 0;