import sys
import getopt
import random
import time
import numpy
# import math    TODO: UNUSED
import arcpy
import csv
from arcpy import env
import os

import my_tools as tools

# Efficiency: 0' 10'' importing packages
#             6' 17'' retrieving day data for January 2010 (500,000 rows read, 1,326 rows/sec), assembly of shapefiles and merging to month file 
#             1' 06'' density analysis, simulating density in 99 random distributions of crime (0.6 secs/simulation) and envelope hypothesis testing

# Note: file names of ArcGIS datasets cannot have more than 13 characters

try:
    opts = getopt.getopt(sys.argv[1:],"m:",["model="])
except getopt.GetoptError:
    print '\nUsage: Exercise3.py -m <true|false>'
    sys.exit(2)
if len(opts[0]) == 0:
    print '\nWarning: \'model\' argument could not be solved. Empirical density will be computed but not modeled'
    modeling_process = False
elif opts[0][0][1] in ['true','True','TRUE']:
    print '\nWarning: \'model\' argument set to True. Empirical density will be computed and modeled'
    modeling_process = True
elif opts[0][0][1] in ['false','False','FALSE']:
    print '\nWarning: \'model\' argument set to False. Empirical density will be computed but not modeled'
    modeling_process = False
else:
    print '\nWarning: \'model\' argument  could not be solved. Empirical density will be computed but not modeled'
    modeling_process = False

arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True

# This sets workspace based on if a specific path exists on computer
if os.path.exists("C:/z_Workspace"):
    print("Setting workspace for home computer.")
    # Sets workspace for home computer
    env.workspace = "C:/z_Workspace/Python/Chicago"
    env.scratchWorkspace = "C:/z_Workspace/Python/scratch"
else:
    print("Setting workspace for work computer.")
    # Sets workspace for work computer
    env.workspace = "C:/Python27/ArcGIS10.4/data/chi2"
    env.scratchWorkspace = "C:/Python27/ArcGIS10.4/data"

# workspace and scratchWorkspace commented out because I have a custom
# workspace settings coded above
'''
env.workspace = "C:\Python27\ArcGIS10.4\data\chi2"  # Set workspace
env.scratchWorkspace = "C:\Python27\ArcGIS10.4\data"
'''

raw_data = "Chicago_crime_counts_2009-2011.csv"
root_name = "chicago_"

'''
Import raw 2010 data to shapefile
'''

year = 2010
segment = 5
progress = 1
sys.stdout.write("\nProcessing days of January, " + str(year) + ": ")
for year_day in range(1,32):  # Only January 2010 is assessed
    percentile = int(year_day/float(31)*100)
    if percentile/segment >= progress:
        sys.stdout.write(str(percentile/segment*segment) + '% ')
        progress = percentile/segment + 1
    root_yday_name = root_name + "D" + str(year_day)
    with open(env.workspace + "\\" + root_yday_name + ".csv", 'w') as out_csv:
        writer = csv.writer(out_csv, delimiter=',',quotechar='"')
        writer.writerow(["ID","Case","Date","Type","X","Y","Lat","Long"])
        with open(env.workspace + "\\" + raw_data, 'r') as raw_csv:
            reader = csv.reader(raw_csv, delimiter=',', quotechar='"')
            next(raw_csv)
            for row in reader:
                time_structure = time.strptime(row[2],"%m/%d/%Y %I:%M:%S %p") ## TO-DO: Print the object and study its data
                if time_structure[0] == year and time_structure[7] == year_day and row[3] == 'NARCOTICS':
                    writer.writerow(row)
    arcpy.MakeXYEventLayer_management(env.workspace + "\\" + root_yday_name + ".csv", "X", "Y", root_yday_name, "", "")
    arcpy.FeatureClassToFeatureClass_conversion(root_yday_name, env.workspace, root_yday_name)

'''
Merge weekly crime data by month
'''

map_names = arcpy.ListFeatureClasses(root_name + "D*")
month_length = [31,28,31,30,31,30,31,31,30,31,30,31]
yday = 1
month = 1
while yday < 32: # Only January 2010 is assessed
    maps_to_merge = []
    indexes_to_use = range(yday-1,yday-1+month_length[month-1])
    for i in indexes_to_use:
        maps_to_merge.append(map_names[i])
    arcpy.Merge_management(maps_to_merge,root_name + "M" + str(month))
    yday += month_length[month-1]
    month += 1

'''
Compute monthly density and significance of local density singals with Monte Carlo simulation
'''

foot = 0.3048
cell_size = 300/foot                                                              # Side cell resolution
random.seed(time.time())                    
chicago_neighborhoods = 'nhoods/Neighborhoods_2012b.shp'                                 # Polygon shapefile with representing the extent of the study area
city_extent = arcpy.Dissolve_management(chicago_neighborhoods,"chicago_extent")
city_extent_rast = arcpy.FeatureToRaster_conversion(city_extent, "Id", "chicago_ext", cell_size)
n_simulations = 199                                                               # Ideally 199 at least for two-tailed alpha = 0.01
alpha = 0.01                                                                      # Portion in the envelope of simulated intensities to be deemed anomalous, as alpha == percent/100
mbr_poly = arcpy.MinimumBoundingGeometry_management(city_extent, "chicago_mbr", "ENVELOPE", "ALL")
mbr_rast = arcpy.FeatureToRaster_conversion(mbr_poly, "Id", "chicago_mbr", cell_size)
min_x = float(arcpy.GetRasterProperties_management(mbr_rast,"LEFT").getOutput(0))
max_x = float(arcpy.GetRasterProperties_management(mbr_rast,"RIGHT").getOutput(0))
min_y = float(arcpy.GetRasterProperties_management(mbr_rast,"BOTTOM").getOutput(0))
max_y = float(arcpy.GetRasterProperties_management(mbr_rast,"TOP").getOutput(0))
n_cols = int(arcpy.GetRasterProperties_management(mbr_rast,"COLUMNCOUNT").getOutput(0))
n_rows = int(arcpy.GetRasterProperties_management(mbr_rast,"ROWCOUNT").getOutput(0))
n_cells = n_cols * n_rows
envelope = numpy.zeros((n_simulations+1,n_cells))
        
for month in range(1,2):  # Only January 2010 is assessed
        
	'''
        Compute local point density with 2-dimensional isotropic Gaussian  kernel of 1000-meter radius
        '''

        input_data = root_name + "M" + str(month) + ".shp"
        n_monthly_events = arcpy.GetCount_management(input_data)
        population_field = "NONE"
        output_intensity = arcpy.sa.KernelDensity(env.workspace + "\\" + input_data, population_field, cell_size,1000/foot,'SQUARE_FEET')
        output_map = "lambda_M" + str(month)  
        output_intensity.save(env.workspace + "\\" + output_map)

        if modeling_process:
    	    '''
            Test anomaly of empirical local density compared to random distributions of crime events
            '''
            sys.stdout.write("\n\nProcessing simulations: ")
            segment = 5
            progress = 1
	    for simulation in range(1,n_simulations+1):
                percentile = int(simulation/float(n_simulations)*100)
                if percentile/segment >= progress:
                    sys.stdout.write(str(percentile/segment*segment) + '% ')
                    progress = percentile/segment + 1
	        arcpy.CreateRandomPoints_management(env.workspace,"simulation",city_extent,"",n_monthly_events)
	        sim_intensity = arcpy.sa.KernelDensity(env.workspace + "\\" + 'simulation.shp', population_field, cell_size, 1000/foot, 'SQUARE_FEET')
                envelope[simulation] = arcpy.RasterToNumPyArray(sim_intensity,arcpy.Point(min_x, min_y),n_cols,n_rows).ravel()  # Append simulated intensity to the envelope of posible simulations
            envelope[-1] = arcpy.RasterToNumPyArray(output_intensity,arcpy.Point(min_x, min_y),n_cols,n_rows).ravel()           # Append empirical intensity to the envelope of posible simulations

            pseudoalpha = numpy.zeros(n_cells)                                                              # For each cell in output map:
            CSR_test = numpy.zeros(n_cells,dtype='int')
            for cell in range(n_cells):
                pseudoalpha[cell] = 1-sum(envelope[:,cell] < envelope[-1, cell])/float(n_simulations+1)     # Compute the rank of the empirical intensity with regard to simulated intensities
                if pseudoalpha[cell] <= float(alpha)/2 or pseudoalpha[cell] >= 1-float(alpha)/2:            # Solve the test based on the 1-alpha confidence interval
                    CSR_test[cell] = 1

            pseudoalpha = numpy.split(pseudoalpha,n_rows)
            pseudoalpha = numpy.array(pseudoalpha)
            pseudoalpha_raster = arcpy.NumPyArrayToRaster(pseudoalpha,
                                                          lower_left_corner=arcpy.Point(min_x, min_y),
                                                          x_cell_size=cell_size)
            pseudoalpha_raster = pseudoalpha_raster + city_extent_rast
            pseudoalpha_raster.save("pseudoalpha")
        
            CSR_test = numpy.split(CSR_test,n_rows)
            CSR_test = numpy.array(CSR_test)
            CSR_test_raster = arcpy.NumPyArrayToRaster(CSR_test,
                                                       lower_left_corner=arcpy.Point(min_x, min_y),
                                                       x_cell_size=cell_size)
            CSR_test_raster = CSR_test_raster + city_extent_rast
            CSR_test_raster.save("CSR_test")
            raster_table = arcpy.BuildRasterAttributeTable_management(CSR_test_raster, "Overwrite")
            n_cells = 0
            with arcpy.da.SearchCursor(raster_table, ['VALUE', 'COUNT']) as cursor:
                for row in cursor:
                    n_cells += row[1]
                    if row[0] == 1:
                        alternative_space = row[1]
            print '\n\nFrequency of anomalous space (alpha = ' + str(alpha) + ' | simulations = ' + str(n_simulations) + '): ' + str(round(alternative_space/float(n_cells),4))

print '\nDone'