function out=hedgeGeography(mu,va,weather,stationData,fitModels,n)
% hedgeGeography reads in geographic weather data, picks a location at
% random, attempts to fit b and d parameters for the seasonal weather at
% that location. If the parameters are fit, runs a bet-hedging scenario at
% this location and collects summary data of this location and the model
% performance here.
%=========================================================================
% Usage:
% out=hedgeGeography(mu,va,weather,stationData,fitModels,n)
% Inputs:
%       * mu = mean light choice probability
%       * va = variance in light choice probability
%       * weather = struct of weather data, e.g. wp in hedgeWeatherPack.mat
%       * stationData = structure containing weather station locations and
%            seasonal temperature normal data. Variable "stationData" in
%            hedgeStationData.mat
%       * fitModels = cell array of previously fit station data (i.e. prior
%            output of this function)
%       * n = number of stations to attempt to fit
% Output, a cell array, with rows = stations and columns as follows
%       1) output of hedgeAnalytic for the calibrated adaptive tracking strategy
%       2) output of hedgeAnalytic for the calibrated bet-hedging strategy
%       3) final population size of the bet-hedging strategy
%       4) latitude of station
%       5) longitude of station
%       6) season length in days
%       7) mean seasonal temperature
%       8) max temperature differential across season
%       9) calibrated b parameter
%       10)  calibrated d parameter

%unpack list of stations, shuffle it, take first n as working list
stationList=stationData.stationLatLongLabels;
stationPerm=randperm(length(stationList));
stationList=stationList(stationPerm);
stationList=stationList(1:n);

%initialize output
out=cell(n,10);

%pull out weather data of stations already (array rows 6-8), and their
%corresponding b and d parameters
knownSeasons=cell2mat(fitModels(:,6:10));

% for number of attempted stations
for i=1:n
    
    %get name of current station, delete second character for NOAA
    %comparability
    targetStation=stationList{i};
    targetStation(2)=[];
    
    %scan for this station among list of all stations, break when found
    for j=1:length(stationData.stationTemperatureLabels)
        if isequal( stationData.stationTemperatureLabels{j},targetStation)
            break;
        end
    end
    
    %build the annual temp data month by month
    normals=[];
    normals=[normals stationData.stationTemperatureData(j,2:32)];
    normals=[normals stationData.stationTemperatureData(j+1,2:29)];
    normals=[normals stationData.stationTemperatureData(j+2,2:32)];
    normals=[normals stationData.stationTemperatureData(j+3,2:31)];
    normals=[normals stationData.stationTemperatureData(j+4,2:32)];
    normals=[normals stationData.stationTemperatureData(j+5,2:31)];
    normals=[normals stationData.stationTemperatureData(j+6,2:32)];
    normals=[normals stationData.stationTemperatureData(j+7,2:32)];
    normals=[normals stationData.stationTemperatureData(j+8,2:31)];
    normals=[normals stationData.stationTemperatureData(j+9,2:32)];
    normals=[normals stationData.stationTemperatureData(j+10,2:31)];
    normals=[normals stationData.stationTemperatureData(j+11,2:32)];
    % convert from °Fx10 to °C
    normals=(normals/10-32)*5/9;
    
    %if it gets to be more than 10°C at least once
    if sum(normals>=10)>1
        seasonStart=find(normals>=6.5);     %find first day of >6.5°C
        seasonStart=seasonStart(1);
        seasonEnd=find(normals>=10);        %find last day of <10°C
        seasonEnd=seasonEnd(end);
    else
        continue;
    end
    
    %define season by those endpoints
    normals=normals(seasonStart:seasonEnd);
    
    %stats on this season
    normalsStats=[length(normals) mean(normals) max(normals)-min(normals)];
    
    %position of this station
    lati=stationData.stationLatLongs(stationPerm(i),1);
    longi=stationData.stationLatLongs(stationPerm(i),2);
    
    %check if this station is already calibrated (start by assuming no)
    knownBool=0;
    stationStats=[lati longi normalsStats];         %5 values for comparison
    knownStats=cell2mat(fitModels(:,4:8));          %list of corresponding known values
    %see if any of the known vectors match that at hand
    for j=1:size(knownStats,1); if isequal(stationStats,knownStats(j,:)); knownBool=1; end; end;
    %if a match, indicate this, and go to next station
    if knownBool==1; disp('Already Solved!'); continue; end;
    
    % find the closest season in the known data base
    seasonZData=[knownSeasons(:,1:3);normalsStats];     %known season values with season at hand season values tacked to the end
    seasonZData=zscore(seasonZData);                    %normalize by column
    seasonSearch=seasonZData(end,:);                    %search vector is now-normalized last row
    seasonZData(end,:)=[];                              %strip it off known matrix
    seasonMatches=zeros(size(seasonZData,1),1);         %initialize search score output
    %get euclidean distance from this season's values to known seasons
    for j=1:size(seasonZData,1); seasonMatches(j)=sqrt(sum((seasonSearch-seasonZData(j,:)).^2)); end   
    %find the closest match season
    whichIsBestMatch=find(seasonMatches==min(seasonMatches));
    whichIsBestMatch=whichIsBestMatch(1);
    
    %display the iteration number and stats for this season and closest
    %known season. Comment these two lines in if this output is desired
%     disp(i);                                                                                                  
%     disp([ length(normals) mean(normals) max(normals)-min(normals); knownSeasons(whichIsBestMatch,1:3)]);
    
    %calibrate the b and d parameters under a adaptive tracking assumption
    calibModel=hedgeBDCalibrate(mu,va,normals,weather,knownSeasons(whichIsBestMatch,4),knownSeasons(whichIsBestMatch,5));
    
    %extract the weather data
    hedgeWeatherPack.tempHist=normals;
    hedgeWeatherPack.cloudCover=normals./normals*0;
    
    %run the model under the bet-hedging strategy, and the extracted
    %weather data
    hedgeRun=hedgeAnalytic(0,weather,mu,va,calibModel.b,calibModel.d,0,0,7,0,'pre',hedgeWeatherPack);
    
    %if the model attained a very low value (and perhaps a numerically
    %unstable value)
    if min(sum(hedgeRun.pops))<.001
        for j=1:10; out{i,j}=NaN; end;  %ignore this results
    else %otherwise output the pieces into the final cell array
        out{i,1}=calibModel;
        out{i,2}=hedgeRun;
        out{i,3}=sum(hedgeRun.pops(:,end))-100;
        out{i,4}=lati;
        out{i,5}=longi;
        out{i,6}=length(calibModel.normalSeason);
        out{i,7}=mean(calibModel.normalSeason);
        out{i,8}=max(calibModel.normalSeason)-min(calibModel.normalSeason);
        out{i,9}=calibModel.b;
        out{i,10}=calibModel.d;
    end
    
    
end