#!/usr/bin/env python3
"""Calculate 22 base paths and 54 principal sensitivity runs using full grid paths.

Requires Python, NumPy, Pillow and a C++17 compiler. The search preserves the
terrain model's edge costs, graph order, binary64 arithmetic and heap tie break.
The 56 shared baseline scenarios are checked against the retained model costs
and display polylines. The C++ helper accelerates the Dijkstra loop; the Python
terrain module builds the grids, calculates edge costs and expands grid steps.
Run from any directory; --original can select another model dependency folder.
"""
from pathlib import Path
import argparse,csv,hashlib,importlib.util,json,math,platform,subprocess,sys,time
import numpy as np

HERE=Path(__file__).resolve().parent
DEFAULT=HERE.parent/'terrain_model'
p=argparse.ArgumentParser();p.add_argument('--original',type=Path,default=DEFAULT);p.add_argument('--out',type=Path,default=HERE);args=p.parse_args()
ORIGINAL=args.original.resolve();OUT=args.out.resolve();OUT.mkdir(parents=True,exist_ok=True)
def log(s):print(time.strftime('%H:%M:%S'),s,flush=True)
def dump(path,data):path.write_text(json.dumps(data,ensure_ascii=False,indent=2)+'\n')
def sha(p):return hashlib.sha256(p.read_bytes()).hexdigest()
sys.dont_write_bytecode=True
script=ORIGINAL/'generate_least_cost_paths.py'
spec=importlib.util.spec_from_file_location('original_model',script);m=importlib.util.module_from_spec(spec);sys.modules[spec.name]=m;spec.loader.exec_module(m)
original_fetch=m.fetch_tile
def fetch_only_cached(cache,zoom,x,y):
 tile=cache/f'z{zoom}_{x}_{y}.png'
 if not tile.is_file():raise FileNotFoundError(f'Archived tile missing; no network request permitted: {tile}')
 return original_fetch(cache,zoom,x,y)
m.fetch_tile=fetch_only_cached
inputs=json.loads((ORIGINAL/'least_cost_path_inputs.json').read_text())
primary=json.loads((ORIGINAL/'least_cost_path_outputs.json').read_text())
sensitivity=json.loads((ORIGINAL/'least_cost_path_sensitivity.json').read_text())
point_index={p['key']:m.Point(p['key'],p['name'],*p['coords']) for p in inputs['points']}
configs=[dict(key='primary',zoom=8,downsample=4,margin_deg=.9,neighbourhood_radius=2,terrain_weight_scale=1.)]+list(m.SENSITIVITY_VARIANTS)
archived=[primary]+sensitivity['runs']
principal_keys=[f'{site}-{source}' for site in ('harb','karahan','gobekli') for source in ('central','bingol','nemrut')]
original_routes={r['key']:r for r in primary['routes']}
matrix_routes=[]
for site,target in [('harb','harbetsuvan'),('karahan','karahan'),('gobekli','gobekli'),('sefer','sefer'),('tasli','tasli'),('hamzan','hamzan')]:
 for suffix,source in [('central','gollu'),('bingol','bingol'),('nemrut','nemrut')]:
  routekey=f'{site}-{suffix}'
  if site in ('sefer','hamzan') and suffix=='bingol':routekey=f'{site}-east'
  matrix_routes.append(dict(key=routekey,source_key=source,target_key=target,analysis_role='principal' if site in ('harb','karahan','gobekli') else 'supplementary',endpoint_status='illustrative' if site in ('tasli','hamzan') else 'original model endpoint'))
matrix_routes.extend({**{k:r[k] for k in ('key','source_key','target_key')},'analysis_role':'regional comparator','endpoint_status':'original model endpoint'} for r in primary['routes'] if r['target_key'] in ('jerf','kortik'))
assert len(matrix_routes)==22 and len({r['key'] for r in matrix_routes})==22
dump(OUT/'complete_matrix_inputs.json',{'points':inputs['points'],'routes':matrix_routes,'configurations':configs,'scope':'Three selected source-zone alternatives at each of six unchanged local endpoints, plus four unchanged regional comparator paths. Not an exhaustive test of possible sources.','original_input_sha256':sha(ORIGINAL/'least_cost_path_inputs.json')})
fields=['elev','lats','lons','east_west','north_south','ruggedness','topographic_position','flow_accumulation']
cache=OUT/'calculation_cache';cache.mkdir(exist_ok=True)
exe=HERE/'exact_dijkstra'
if not exe.exists() or exe.stat().st_mtime < (HERE/'exact_dijkstra.cpp').stat().st_mtime:
 subprocess.run(['clang++','-std=c++17','-O2','-fno-fast-math','-ffp-contract=off',str(HERE/'exact_dijkstra.cpp'),'-o',str(exe)],check=True)
offsets=[(r,c) for r in (-1,0,1) for c in (-1,0,1) if r or c]
results=[];metrics=[];manifest={str(p.relative_to(ORIGINAL)):sha(p) for p in sorted(ORIGINAL.rglob('*')) if p.is_file() and '__pycache__' not in p.parts and p.suffix!='.pyc'}
fingerprint=cache/'original_fingerprint.json'
if fingerprint.exists():
 assert json.loads(fingerprint.read_text())==manifest, 'Calculation cache belongs to different original files; use a fresh --out directory.'
else:dump(fingerprint,manifest)
assert len(configs)==len(archived)==7
for cfg,archive in zip(configs,archived):
 key=cfg['key']
 if key!='primary':assert archive['run_label']==key
 route_specs=matrix_routes if key=='primary' else [r for r in matrix_routes if r['key'] in principal_keys]
 old_index={r['key']:r for r in archive['routes']}
 log(f'Start {key} ({len(route_specs)} scenarios)')
 used={k for route in archive['routes'] for k in (route['source_key'],route['target_key'])}
 model_points=[point_index[p['key']] for p in inputs['points'] if p['key'] in used]
 xs,ys=m.tile_window(model_points,cfg['zoom'],cfg['margin_deg'])
 signature=f'z{cfg["zoom"]}_f{cfg["downsample"]}_x{xs[0]}_{xs[-1]}_y{ys[0]}_{ys[-1]}'
 gridfile=cache/f'{signature}.npz'
 if gridfile.exists():
  z=np.load(gridfile);grid=m.CostGrid(**{f:z[f] for f in fields},tile_ids=z['tile_ids'].tolist())
 else:
  grid=m.build_cost_grid(model_points,cfg['zoom'],cfg['downsample'],cfg['margin_deg'],ORIGINAL/'terrain_cache')
  np.savez_compressed(gridfile,**{f:getattr(grid,f) for f in fields},tile_ids=np.array(grid.tile_ids))
 rows,cols=grid.elev.shape;log(f'Grid {rows} x {cols}; build/reuse exact Python edge costs')
 edgefile=cache/f'{signature}_weight{cfg["terrain_weight_scale"]}.f64'
 if not edgefile.exists():
  edges=np.full((rows,cols,8),np.inf,dtype=np.float64)
  for row in range(rows):
   for col in range(cols):
    for i,(dr,dc) in enumerate(offsets):
     nr,nc=row+dr,col+dc
     if 0<=nr<rows and 0<=nc<cols:edges[row,col,i]=m.edge_cost(grid,row,col,nr,nc,cfg['terrain_weight_scale'])
  edges.tofile(edgefile)
  del edges
 queries=[]
 for route in route_specs:
  src=point_index[route['source_key']];dst=point_index[route['target_key']]
  start=m.nearest_cell(src.lat,src.lon,grid.lats,grid.lons);goal=m.nearest_cell(dst.lat,dst.lon,grid.lats,grid.lons)
  queries.append((start[0]*cols+start[1],goal[0]*cols+goal[1]))
 queryfile=cache/f'{key}_queries.txt';queryfile.write_text(''.join(f'{a} {b}\n' for a,b in queries))
 rawfile=cache/f'{key}_search.json'
 if not rawfile.exists():subprocess.run([str(exe),str(rows),str(cols),str(cfg['neighbourhood_radius']),str(edgefile),str(queryfile),str(rawfile)],check=True)
 raw=json.loads(rawfile.read_text());assert len(raw)==len(route_specs);run={'configuration':cfg,'grid_shape':[rows,cols],'grid_extent':archive['grid_extent'],'tile_ids':grid.tile_ids,'routes':[]}
 for route,found in zip(route_specs,raw):
  routekey=route['key'];source=point_index[route['source_key']];target=point_index[route['target_key']]
  old=old_index.get(routekey)
  parent=[divmod(i,cols) for i in found['indices']]
  expanded=[parent[0]]
  segment_group_costs=[]
  for a,b in zip(parent,parent[1:]):
   segments=m.bresenham_cells(*a,*b);expanded.extend(segments[1:])
   segment_group_costs.append(sum(m.edge_cost(grid,*x,*y,cfg['terrain_weight_scale']) for x,y in zip(segments,segments[1:])))
  recomputed_cost=sum(segment_group_costs)
  pc=[(float(grid.lats[r]),float(grid.lons[c])) for r,c in parent]
  ec=[(float(grid.lats[r]),float(grid.lons[c])) for r,c in expanded]
  display=m.force_exact_endpoints(m.path_to_latlon(parent,grid.lats,grid.lons),source,target)
  src_connector=m.haversine_km(source.lat,source.lon,*ec[0]);dst_connector=m.haversine_km(*ec[-1],target.lat,target.lon)
  gridlen=0.
  for (r,c),(nr,nc) in zip(expanded,expanded[1:]):
   if r==nr:d=float(grid.east_west[r])
   elif c==nc:d=float(grid.north_south[min(r,nr)])
   else:d=math.hypot(float(grid.east_west[min(r,nr)]),float(grid.north_south[min(r,nr)]))
   gridlen+=d
  total=gridlen+src_connector+dst_connector
  display_match=[[a,b] for a,b in display]==old['path'] if old else None
  rec={
   'run':key,'key':routekey,'source':source.name,'target':target.name,'source_key':source.key,'target_key':target.key,
   'analysis_role':route['analysis_role'],'endpoint_status':route['endpoint_status'],'historical_scenario':old is not None,
   'parent_vertex_count':len(parent),'expanded_vertex_count':len(expanded),
   'grid_path_length_km':gridlen,'source_connector_km':src_connector,'target_connector_km':dst_connector,'total_path_length_km':total,
   'expanded_haversine_grid_length_km':m.polyline_length(ec),'expanded_haversine_total_length_km':m.polyline_length(ec)+src_connector+dst_connector,
   'parent_chain_haversine_grid_length_km':m.polyline_length(pc),'parent_chain_total_length_km':m.polyline_length(pc)+src_connector+dst_connector,
   'archived_display_length_km':old['length_km'] if old else None,'reproduced_display_length_km':m.polyline_length(display),
   'length_increase_km':total-old['length_km'] if old else None,'length_increase_pct':100*(total/old['length_km']-1) if old else None,
   'archived_composite_cost':old['composite_cost_index'] if old else None,'reproduced_composite_cost':found['cost'],'costed_edge_sum':recomputed_cost,
   'cost_difference':found['cost']-old['composite_cost_index'] if old else None,'display_geometry_exact_match':display_match,
   'display_length_difference':m.polyline_length(display)-old['length_km'] if old else None,
   'straight_line_km':m.haversine_km(source.lat,source.lon,target.lat,target.lon),
  }
  rec['excess_over_straight_line_pct']=100*(total/rec['straight_line_km']-1)
  metrics.append(rec)
  full={**rec,'parent_grid_cells':parent,'expanded_grid_cells':expanded,'expanded_grid_coordinates':ec,'source_coordinates':[source.lat,source.lon],'target_coordinates':[target.lat,target.lon],'archived_display_coordinates':display}
  run['routes'].append(full)
  log(f'{routekey}: {total:.6f} km; historical={old is not None}; display exact={display_match}; cost={found["cost"]:.9f}')
  if (old and (not display_match or abs(rec['cost_difference'])>1e-10 or abs(rec['display_length_difference'])>1e-10)) or abs(found['cost']-recomputed_cost)>1e-10:
   dump(OUT/f'{key}_REPRODUCTION_FAILURE.json',run);raise AssertionError(f'Original results not reproduced: {key}/{routekey}')
 dump(OUT/f'{key}_complete_paths.json',run);results.append(run)
 dump(OUT/'distance_metrics_checkpoint.json',metrics)

def writecsv(path,records):
 with path.open('w',newline='') as f:
  w=csv.DictWriter(f,fieldnames=list(records[0]));w.writeheader();w.writerows(records)
assert len(metrics)==76 and sum(r['historical_scenario'] for r in metrics)==56
writecsv(OUT/'distance_metrics_all_76.csv',metrics)
dump(OUT/'distance_metrics_all_76.json',metrics)
base={r['key']:r for r in metrics if r['run']=='primary'}
contrasts=[]
for cfg in configs:
 ri={r['key']:r for r in metrics if r['run']==cfg['key']}
 for near,far in [(f'{site}-{source}',f'{site}-central') for site in ('harb','karahan','gobekli') for source in ('bingol','nemrut')]:
  a,b=ri[near],ri[far]
  contrasts.append({'run':cfg['key'],'near_route':near,'reference_route':far,'new_total_reduction_pct':100*(1-a['total_path_length_km']/b['total_path_length_km']),'grid_only_reduction_pct':100*(1-a['grid_path_length_km']/b['grid_path_length_km']),'old_display_reduction_pct':100*(1-a['archived_display_length_km']/b['archived_display_length_km']) if a['historical_scenario'] else None})
writecsv(OUT/'source_contrasts.csv',contrasts);dump(OUT/'source_contrasts.json',contrasts)
changes=[]
for r in metrics:
 if r['run']=='primary':continue
 b=base[r['key']]
 changes.append({'run':r['run'],'key':r['key'],'new_total_length_km':r['total_path_length_km'],'baseline_new_total_length_km':b['total_path_length_km'],'relative_change_pct':100*(r['total_path_length_km']/b['total_path_length_km']-1),'cost_change':r['reproduced_composite_cost']-b['reproduced_composite_cost']})
writecsv(OUT/'sensitivity_length_changes.csv',changes);dump(OUT/'sensitivity_length_changes.json',changes)
summary={
 'number_of_calculated_scenarios':len(metrics),'number_of_base_scenarios':22,'number_of_principal_scenarios':9,'number_of_sensitivity_scenarios':54,'number_of_reproduced_scenarios':56,'all_historical_display_geometry_matches':all(r['display_geometry_exact_match'] for r in metrics if r['historical_scenario']),'maximum_cost_difference':max(abs(r['cost_difference']) for r in metrics if r['historical_scenario']),'maximum_display_length_difference':max(abs(r['display_length_difference']) for r in metrics if r['historical_scenario']),
 'primary_principal_rows':[base[k] for k in principal_keys],
 'primary_contrasts':[r for r in contrasts if r['run']=='primary'],
 'sensitivity_total_length_change_pct_min':min(r['relative_change_pct'] for r in changes),'sensitivity_total_length_change_pct_max':max(r['relative_change_pct'] for r in changes),
 'sensitivity_ordering_retained':all(r['new_total_reduction_pct']>0 for r in contrasts),
 'three_by_three_total_changes':[r for r in changes if r['run']=='three_by_three_neighbourhood'],
 'original_files_unchanged':all(sha(ORIGINAL/p)==h for p,h in manifest.items()),
 'source_file_sha256':manifest,'runtime':{'python':platform.python_version(),'numpy':np.__version__},
}
dump(OUT/'analysis_summary.json',summary)
dump(OUT/'distance_metric_definitions.json',{
 'scope':'Complete three-source by six-local-endpoint matrix plus four regional comparator paths, preserving original grids, endpoints, cost assumptions and the corrected full-path distance metric. Nine principal paths are tested in each of six sensitivity variants. All 56 historical calculations are reproduced.',
 'grid_path_length_km':'Sum of the horizontal distance of every elementary step actually costed, after expanding each parent link with the original Bresenham function. Cardinal distances use the original float32 row distance tables; diagonal distances use hypot(original east-west, original north-south) at the lesser row index, exactly as original edge_cost.',
 'source_connector_km':'Haversine distance from exact source point to source grid-cell centre; always included, even below the original display insertion threshold of 0.01 km.',
 'target_connector_km':'Haversine distance from target grid-cell centre to exact target point; always included, even below 0.01 km.',
 'total_path_length_km':'grid_path_length_km + source_connector_km + target_connector_km. Recommended reported distance. Connectors are geometric and were not costed or used in path optimisation.',
 'expanded_haversine_grid_length_km':'Diagnostic sum of haversine distances along expanded elementary coordinates; distinguishes great-circle measurement from the diagonal grid-distance approximation used in original costs.',
 'parent_chain_total_length_km':'Diagnostic unsimplified parent-chain haversine length plus both connectors. This still omits intermediate cells in 5×5 moves and is not the full costed grid path.',
 'archived_display_length_km':'Historical sum of haversine lengths of simplified display line and original thresholded endpoint handling.',
 'percentage_rule':'100*(1-nearer total_path_length_km/reference total_path_length_km), using unrounded binary64 results.',
 'reproduction':'All original-Python edge costs are precomputed without altered arithmetic. Compiled Dijkstra retains the original neighbour order, lexicographic (score,index) queue order, closed-set logic, strict relaxation and binary64 additions. No fast-math or floating-point reassociation. All 56 archived display polylines, lengths and accumulated costs must reproduce before output is accepted.'})
log('DONE: 22 base and 54 sensitivity scenarios calculated; all 56 historical results reproduced.')