# -*- coding: utf-8 -*-
"""
Created on Thu Dec  6 15:33:49 2018

@author: Nikolina Mileva
"""

import zarr
import numpy as np
import dask.array as da
from dask.diagnostics import ProgressBar
from parameters import (windowSize, logWeight, temp, mid_idx, numberClass, spatImp, 
                        specUncertainty, tempUncertainty, path)

 


# Flatten blocks inside a dask array            
def block2row(array, row, folder, block_id=None):
    if array.shape[0] == windowSize:
        # Parameters	
        name_string = str(block_id[0] + 1)
        m,n = array.shape
        u = m + 1 - windowSize
        v = n + 1 - windowSize

    	# Get Starting block indices
        start_idx = np.arange(u)[:,None]*n + np.arange(v)

    	# Get offsetted indices across the height and width of input array
        offset_idx = np.arange(windowSize)[:,None]*n + np.arange(windowSize)

    	# Get all actual indices & index into input array for final output
        flat_array = np.take(array,start_idx.ravel()[:,None] + offset_idx.ravel())

        # Save to (dask) array in .zarr format
        file_name = path + folder + name_string + 'r' + row + '.zarr'
        zarr.save(file_name, flat_array)
    
    return array


# Divide an image in overlapping blocks   
def partition(image, folder):
    image_da = da.from_array(image, chunks = (windowSize,image.shape[1]))
    image_pad = da.pad(image_da, windowSize//2, mode='constant')
    
    for i in range(0,windowSize):
        row = str(i)
        block_i = image_pad[i:,:]
        block_i_da = da.rechunk(block_i, chunks=(windowSize,image_pad.shape[1]))
        block_i_da.map_blocks(block2row, dtype=int, row=row, folder=folder).compute()


# Create a list of all files in the folder and stack them into one dask array
def da_stack(folder, shape):
    da_list = [] 
    full_path = path + folder
    max_blocks = shape[0]//windowSize + 1 
    
    for block in range(1,max_blocks + 1):
        for row in range(0,windowSize):
            name = str(block) + 'r' + str(row)
            full_name = full_path + name + '.zarr'
            try:
                da_array = da.from_zarr(full_name)
                da_list.append(da_array) 
            except Exception:
                continue
      
    return da.rechunk(da.concatenate(da_list, axis=0), chunks = (shape[1],windowSize**2))


# Calculate the spectral distance
def spectral_distance(fine_image_t0, coarse_image_t0):
    spec_diff = fine_image_t0 - coarse_image_t0
    spec_dist = 1/(abs(spec_diff) + 1.0)
    print ("Done spectral distance!", spec_dist)
    
    return spec_diff, spec_dist


# Calculate the temporal distance    
def temporal_distance(coarse_image_t0, coarse_image_t1):
    temp_diff = coarse_image_t1 - coarse_image_t0
    temp_dist = 1/(abs(temp_diff) + 1.0)
    print ("Done temporal distance!", temp_dist)
    
    return temp_diff, temp_dist
   

# Calculate the spatial distance    
def spatial_distance(array):
    coord = np.sqrt((np.mgrid[0:windowSize,0:windowSize]-windowSize//2)**2)
    spat_dist = np.sqrt(((0-coord[0])**2+(0-coord[1])**2))
    rel_spat_dist = spat_dist/spatImp + 1.0 # relative spatial distance
    rev_spat_dist = 1/rel_spat_dist # relative spatial distance reversed
    flat_spat_dist = np.ravel(rev_spat_dist)
    spat_dist_da = da.from_array(flat_spat_dist, chunks=flat_spat_dist.shape)
    print ("Done spatial distance!", spat_dist_da)
    
    return spat_dist_da


# Define the threshold used in the dynamic classification process
def similarity_threshold(fine_image_t0):#, st_dev):
    fine_image_t0 = da.where(fine_image_t0==0, np.nan, fine_image_t0)
    st_dev = da.nanstd(fine_image_t0, axis=1)# new
    sim_threshold = st_dev*2/numberClass 
    print ("Done similarity threshold!", sim_threshold)

    return sim_threshold


# Define the spectrally similar pixels within a moving window    
def similarity_pixels(fine_image_t0):
    sim_threshold = similarity_threshold(fine_image_t0)
    # possible to implement as sparse matrix
    similar_pixels = da.where(abs(fine_image_t0 - 
                                  fine_image_t0[:,mid_idx][:,None])
        <= sim_threshold[:,None], 1, 0) #sim_threshold[:,mid_idx][:,None], 1, 0) # new
    print ("Done similarity pixels!", similar_pixels)
   
    return similar_pixels
        

# Apply filtering on similar pixels 
def filtering(fine_image_t0, spec_dist, temp_dist, spec_diff, temp_diff):
    similar_pixels = similarity_pixels(fine_image_t0) 
    max_spec_dist = abs(spec_diff)[:,mid_idx][:,None] + specUncertainty + 1
    max_temp_dist = abs(temp_diff)[:,mid_idx][:,None] + tempUncertainty + 1  
    spec_filter = da.where(spec_dist>1.0/max_spec_dist, 1, 0)
    st_filter = spec_filter
    
    if temp == True:
        temp_filter = da.where(temp_dist>1.0/max_temp_dist, 1, 0)
        st_filter = spec_filter*temp_filter  
        
    similar_pixels_filtered = similar_pixels*st_filter
    print ("Done filtering!", similar_pixels_filtered)

    return similar_pixels_filtered # sim_pixels_sparse
    

# Calculate the combined distance
def comb_distance(spec_dist, temp_dist, spat_dist):
    if logWeight == True:
        spec_dist = da.log(spec_dist + 1)
        temp_dist = da.log(temp_dist + 1)
    
    comb_dist = da.rechunk(spec_dist*temp_dist*spat_dist, 
                           chunks=spec_dist.chunksize)
    print ("Done comb distance!", comb_dist)
    
    return comb_dist
    
        
# Calculate weights
def weighting(spec_dist, temp_dist, comb_dist, similar_pixels_filtered):
    # Assign max weight (1) when the temporal or spectral distance is zero
    zero_spec_dist = da.where(spec_dist[:,mid_idx][:,None] == 1, 1, 0)
    zero_temp_dist = da.where(temp_dist[:,mid_idx][:,None] == 1, 1, 0)
    zero_dist_mid = da.where((zero_spec_dist == 1), 
                             zero_spec_dist, zero_temp_dist)
    shape = da.subtract(spec_dist.shape,(0,1))
    zero_dist = da.zeros(shape, chunks=(spec_dist.shape[0],shape[1]))
    zero_dist = da.insert(zero_dist, [mid_idx], zero_dist_mid, axis=1)
    weights = da.where((da.sum(zero_dist,1)[:,None] == 1), zero_dist, comb_dist)
    
    # Calculate weights only for the filtered spectrally similar pixels
    weights_filt = weights*similar_pixels_filtered
    
    # Normalize weights
    norm_weights = da.rechunk(weights_filt/(da.sum(weights_filt,1)[:,None]), 
                              chunks = spec_dist.chunksize)
    
    print ("Done weighting!", norm_weights)
    
    return norm_weights


# Derive fine resolution reflectance for the day of prediction 
def predict(fine_image_t0, coarse_image_t0, coarse_image_t1, shape):
    spec = spectral_distance(fine_image_t0, coarse_image_t0)
    spec_diff = spec[0]
    spec_dist = spec[1]
    temp = temporal_distance(coarse_image_t0, coarse_image_t1)
    temp_diff = temp[0] 
    temp_dist = temp[1]
    spat_dist = spatial_distance(fine_image_t0)
    comb_dist = comb_distance(spec_dist, temp_dist, spat_dist)
    similar_pixels = filtering(fine_image_t0, spec_dist, temp_dist, spec_diff, 
                               temp_diff)
    weights = weighting(spec_dist, temp_dist, comb_dist, similar_pixels)    
    pred_refl = fine_image_t0 + temp_diff
    weighted_pred_refl = da.sum(pred_refl*weights, axis=1)   
    prediction = da.reshape(weighted_pred_refl, shape)
    print ("Done prediction!")
    
    return prediction
    
 
# Compute the results (converts the dask array to a numpy array)   
def starfm(fine_image_t0, coarse_image_t0, coarse_image_t1, profile, shape):
    print ('Processing...')
    prediction_da = predict(fine_image_t0, coarse_image_t0, coarse_image_t1, shape)
    with ProgressBar():
         prediction = prediction_da.compute()
    
    return prediction



























        
        
        
