# -*- coding: utf-8 -*-
"""
Created on Fri Oct  4 11:58:08 2019

@author: iam
"""

import numpy as np
import numpy.ma as ma
import scipy.stats as stats

from scipy.spatial import cKDTree
from scipy.spatial import Delaunay
from scipy import spatial
from scipy import interpolate

import collections # for defaultdict

from . import f_models
from . import simplex as six
from . import IS_stat
from .IS_stat import IS
from .candybox import CandyBox
from . import sball 



#
# l2 Voronoi cKDTree IS localization
#
 
# Rbf_estimation so far
def L2_Voronoi_cKDTree(sample_box, nis=50000):    
    nsim, nvar = np.shape(sampled_plan_R)
    # Estimation of pf given a list of red and green points (Voronoi)
    
    sampled_plan_R_ma = np.ma.asarray(sampled_plan_R)
    
    # tady bych nedaval Rd 
    #grad = np.abs(gradient([0 for j in range(nvar)]))
    tree = cKDTree(sampled_plan_R)
    

    rbf_func = interpolate.Rbf(sampled_plan_R[:,0],sampled_plan_R[:,1],Z)

    Rbf_estimation = 0 # inicializace
    L2_Voronoi_cKDTree_test = 0
    #failure_points_indexes = np.argwhere(failsi)
    
    points_probabilities = np.empty(nsim) # probabilities of partial (failure) event
    
    # loop over points (need to integrate red regions)
    for i in range(nsim): # loop over failing points only
        
        sampled_plan_R_ma.mask = ma.nomask
        sampled_plan_R_ma[i] = ma.masked
        
        # neosvědčílo se
#==============================================================================
#         delta_Rd_matrix = sampled_plan_Rd_ma - sampled_plan_Rd[i]
#         mindist = [np.min(np.abs(np.where(delta_Rd_matrix[:,j] < 0,delta_Rd_matrix[:,j], f_i[j].std() ))) + np.min(np.abs(np.where(delta_Rd_matrix[:,j] > 0,delta_Rd_matrix[:,j], f_i[j].std() ))) for j in range(nvar)]
#                                          
#         # set the minimum distance as the standard deviation of IS densisty
#         h_i = [stats.norm(sampled_plan_Rd[i,j], mindist[j] ) for j in range(nvar)] #! dosadit  standard deviation pddle chutí
#==============================================================================
        # find distance to the nearest sampling point (from all points)
        mindist = np.min(np.sum(np.square(sampled_plan_R_ma - sampled_plan_R[i]), axis=1))**0.5
                                 
        # set the minimum distance as the standard deviation of IS densisty
        h_i = [stats.norm(sampled_plan_R[i,j], 2*mindist ) for j in range(nvar)] #! dosadit  standard deviation pddle chutí
                                 
        # use IS sampling density with center equal to the current "red" point
                                 
        # select nis = 100 points from IS density and 
        # if the point has its nearest neighbor any red point from the sampled_plan, 
                                 
        h_plan = np.zeros((nis, nvar))
        for j in range(nvar):
            h_plan[:, j] = h_i[j].ppf(np.random.random(nis)) # realizace váhové funkce náhodné veličiny
    
        # Rozptyl corrected IS
        weights_sim = np.prod([f_i[j].pdf(h_plan[:, j]) / h_i[j].pdf(h_plan[:, j]) for j in range(nvar)], axis=0) # [f1/h1, ..., fn/hn]
    
       
    
        dd, ii = tree.query(h_plan)
    
        Vor_mask = np.where(ii==i, True, False)
        
        zm = rbf_func(h_plan[:, 0][Vor_mask], h_plan[:, 1][Vor_mask])        
        
        points_probabilities[i] = np.sum(weights_sim[Vor_mask][zm<0]) / nis
        L2_Voronoi_cKDTree_test += np.sum(weights_sim[Vor_mask]) / nis
    
    
    Rbf_estimation = np.sum(points_probabilities)
    
    return Rbf_estimation, L2_Voronoi_cKDTree_test
    





# cKDTree
"""
# dd  - The distances to the nearest neighbors 
# ii - The locations of the neighbors in self.data
# k - The list of k-th nearest neighbors to return. 
# If k is an integer it is treated as a list of [1, … k] (range(1, k+1)). 
# Note that the counting starts from 1
# p - Which Minkowski p-norm to use. 
# 1 is the sum-of-absolute-values “Manhattan” distance 2 is the usual Euclidean distance 
# infinity is the maximum-coordinate-difference distance 

dd, ii = tree.query(h_plan_model, k=1, p=p_norm)
"""


#
# Voronoi estimation
#    (gradient tu nechcu)
def Voronoi_tesselation(sample_box, model_space='Rn', sampling_space=None, p_norm=1, budget=20000, callback=None):
    """
    Voronoi estimations 
    budget=20000
    """
    global_stats = {'Voronoi_failure_rate':0}

    
    
    nsim = sample_box.nsim
    
        
    # jsou to informace nejen pro callback 
    # nis upresníme dále
    estimation={'method': "Voronoi_tesselation",  'p_norm':p_norm, 'nis':0, 'nsim':nsim}
    estimation['model_space'] = model_space
    estimation['sampling_space'] = sampling_space
    
    if len(sample_box.failure_points) > 0:
        nis = max(round(budget/len(sample_box.failure_points)), 100)
        estimation['nis'] = nis
    
        # vytahneme ze sample_boxu rozdělení
        f = sample_box.sampled_plan
        
        # já vím, že sample box pokážde failure_points přepočítavá
        failure_points = sample_box.failure_points 
        
        # zde provadím rozdělení na prostor, ve kterém vzorkujem
        # a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
        sampled_plan_model = getattr(sample_box, model_space)
        tree = cKDTree(sampled_plan_model)
        
        
        if sampling_space is None:
            sampling_space = model_space
            # sing like sampling
            sampled_plan_sing = sampled_plan_model
            tree_sampling = tree
        else:
            sampled_plan_sing = getattr(sample_box, sampling_space)
            # narychlo, moc to nepomůže, neboť asi po 16 vzorcích počítá brut forsem 
            tree_sampling = cKDTree(sampled_plan_sing,  compact_nodes=False, balanced_tree=False)
            
        # find distance to the nearest sampling point (from all points)
        dd2, ii2 = tree_sampling.query(sampled_plan_sing[failure_points], k=[2], p=p_norm)
        mindist_sing = dd2.flatten()
        #ii2 = ii2.flatten()
        
        
        
        # loop over points (need to integrate red regions)
        for i, i2 in zip(sample_box.failure_points, range(len(mindist_sing))): # loop over failing points only
            # use IS sampling density with center equal to the current "red" point
            # set the minimum distance as the standard deviation of IS densisty
            h_i = [stats.norm(sampled_plan_sing[i,j], 1.5*mindist_sing[i2]) for j in range(sample_box.nvar)] #! dosadit  standard deviation podle chutí
            h = f_models.UnCorD(h_i)
            
                                     
            # select nis = 100 points from IS density 
            # sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
            h_plan = IS(f, h, space_from_h='R', space_to_f=sampling_space, Nsim=nis) 
            
            # součet váh nemá cenu kontrolovat, jedná tam nebude, nevyjde
            
            
            h_plan_model = getattr(h_plan, model_space)
            dd, ii = tree.query(h_plan_model, k=1, p=p_norm)
            
            cell_probability = np.sum(h_plan.w[ii==i]) / nis
            global_stats['Voronoi_failure_rate'] += cell_probability
            
            # kolbek ↓
            if callback is not None:
                cell_stats = {'Voronoi_failure_rate': cell_probability}
                cell_stats['cell_probability'] = cell_probability
                cell_stats['event'] = 'failure'
                callback(estimation=estimation, nodes=h_plan[ii==i], cell_stats=cell_stats, out_nodes=h_plan[ii!=i])

        
    result = {**estimation, **global_stats}
    # try to put it to the BlackBox
    try:
        sample_box.estimations.append(result)
    except: # asi ani nebudu neuspěch hlasit
        pass
        
    return result
    

    
    
    
#
# Voronoi 2_point estimation
# 
def Voronoi_2_point_estimation(sample_box, model_space='Rn', sampling_space=None, p_norm=1, gradient=None, budget=20000, callback=None):
    """
    Voronoi_2_point estimations 
    budget=20000
    """
    
    # tak, na rovinu
    # нет ножек - нет мультиков
#    if gradient is None:
#        return Voronoi_2_point_cKDTree(sample_box, model_space=model_space,sampling_space=sampling_space,\
#                                                            p_norm=p_norm, budget=budget, callback=callback)

#    if gradient is None:
#        gradient = lambda *_, **__: np.ones(sample_box.nvar)

    if callback is None:
        callback = lambda *_, **__: None
    
    nsim = sample_box.nsim
    # ощень, ощень скромненько
    p_base = 0.3 # chcu, aby behem první iterace p_base tečiček jistě dopadlo do buňky
    nis_base = int((sample_box.nvar+1)/p_base)
    nis = int(max(budget/nsim, nis_base / stats.binom.sf(int(sample_box.nvar+1), nis_base, p_base)))
        
    # jsou to informace pro callback 
    estimation={'method': "Voronoi_2_point_estimation",  'p_norm':p_norm, 'nis':nis, 'nsim':nsim}
    estimation['model_space'] = model_space
    estimation['sampling_space'] = sampling_space
    
    
    # vytahneme ze sample_boxu rozdělení
    f = sample_box.sampled_plan
    
    # já vím, že sample box pokážde failsi přepočítavá
    failsi = sample_box.failsi 
    
    PDF = sample_box.pdf(model_space)
    
    # zde provadím rozdělení na prostor, ve kterém vzorkujem
    # a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
    sampled_plan_model = getattr(sample_box, model_space)
    tree = cKDTree(sampled_plan_model)
    
    
    if sampling_space is None:
        sampling_space = model_space
        # sing like sampling
        sampled_plan_sing = sampled_plan_model
        tree_sampling = tree
    else:
        sampled_plan_sing = getattr(sample_box, sampling_space)
        # narychlo, moc to nepomůže, neboť asi po 16 vzorcích počítá brut forsem 
        tree_sampling = cKDTree(sampled_plan_sing,  compact_nodes=False, balanced_tree=False)
        
    # find distance to the nearest sampling point (from all points)
    dd2, ii2 = tree_sampling.query(sampled_plan_sing, k=[2], p=p_norm)
    mindist_sing = dd2.flatten()
    
    # chcu, aby behem první iterace pulka (p_base) tečiček jistě dopadla do buňky
    base_r = sball.Sball(sample_box.nvar).get_r(1 - p_base)
    # np.sum(abs(stats.norm(0, 10).rvs(1000))<sb.get_r(.5)*10)
    
    

    global_stats = collections.defaultdict(int)
    
    
    
    # zde postupně v cyklu prochazíme všemi body (tj. vzorky)
    # a omezujeme se pouse nejbližšími bodíkama
    # tynhlenstím zajišťujeme disjunktnost 
    # a môžeme všechny nasbírané pravděpodobnosti jednoduše sčítat
    for i in range(nsim): # loop over all points, not failing points only
    
        # use IS sampling density with center equal to the current point
        # set the minimum distance as the standard deviation of IS densisty
        # mindist dělíme dvěma - hranice buňky je na půlcesty
        # a dělíme base_r'kem k získání sigmy. Přemejmenším empiricky mně to vychází
        # (u Norm'u vždy zadaváme směrodatnou odchylku)
        h_i = [stats.norm(sampled_plan_sing[i,j], mindist_sing[i]/2/base_r) for j in range(sample_box.nvar)] 
        h = f_models.UnCorD(h_i)
        
                                 
        # select nis points from IS density 
        # sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
        h_plan = IS(f, h, space_from_h='R', space_to_f=sampling_space, Nsim=nis) 
        
        # součet váh nemá cenu kontrolovat, jedná tam nebude, nevyjde
        
        h_plan_model = getattr(h_plan, model_space)
        dd, ii = tree.query(h_plan_model, k=1, p=p_norm)
        
        # nechám s velkým písmenkem
        Vor_mask = ii==i
        h_plan_model_ma = h_plan_model[Vor_mask]
        weights_sim = h_plan.w[Vor_mask]
        # dd1 jsou vzdalenosti tečiček do centra Voroneho buňky
        dd1 = dd[Vor_mask]
        
        trace = 0
        nis_ma = len(weights_sim)
        
        h_plan_sing_ma = getattr(h_plan, sampling_space)[Vor_mask]
        S_bc = np.cov(h_plan_sing_ma, rowvar=False, bias=True, aweights=weights_sim)
        barycenter = np.average(h_plan_sing_ma, axis=0, weights=weights_sim)
        # effective nis
        nis_eff = nis
        #print("start", i)
        
        while trace < np.trace(S_bc):
            # number of nodes inside the cell matters too
            trace = np.trace(S_bc)# * nis_ma
            #print(S_bc)
            #print(nis_ma)
        
            
            
            
            # matika
            w, v = np.linalg.eig(S_bc)
            
            # use IS sampling density with center equal to the cell's barycenter
            # set the minimum distance as the standard deviation of the IS densisty
            # u stats.norm zadáváme směrodatnou odchylku, to je asi správné
            sigmas = np.sqrt(w) / base_r #(sample_box.nvar+2) #! dosadit  standard deviation podle chutí
            h_i = [stats.norm(0, sigmas[j]) for j in range(sample_box.nvar)]
            # rozdělení ve vlastním prostoru
            # select nis = 100 points from IS density 
            h_L = f_models.UnCorD(h_i)(nis)
            
            # здесь уже так легко не отделаемся. Трансформовать кароно.
            h_plan_bc = (v @ h_L.R.T).T
            h_plan_sing = h_plan_bc + barycenter
            
            # sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
            h_plan_part = f.new_sample(h_plan_sing, sampling_space)
            
            
            # jdeme na ty hustoty
            # mně příjde, že je to legalní
            # sice samply podporujou maskovaní, to je ale drahé
            weights_sim_part = h_plan_part.pdf(sampling_space) / h_L.pdf('R') # snad je to správně
            h_plan.add_sample(CandyBox(h_plan_part, w=weights_sim_part))
            
            # vyfiltrujeme vzorky
            h_plan_model_part = getattr(h_plan_part, model_space)
            dd, ii = tree.query(h_plan_model_part, k=1, p=p_norm)
            
            # parta
            Vor_mask_part = ii==i
            weights_sim_part = weights_sim_part[Vor_mask_part]
            
            nis_ma = len(weights_sim_part)
            
            # zajišťovat Vor_mask je docela zbytečně, je to jen pro out_nodes,
            # které se zatím nikdě nepouživá
            Vor_mask = np.append(Vor_mask, Vor_mask_part)
            
            
            
            
            h_plan_model_ma = np.vstack((h_plan_model_ma, h_plan_model_part[Vor_mask_part]))
            weights_sim = np.append(weights_sim, weights_sim_part)
            # dd1 jsou vzdalenosti tečiček do centra Voroneho buňky
            dd1 = np.append(dd1, dd[Vor_mask_part])
            
            # zkusím těžiště počitat jen pro partu - možná tak algoritmus bude agresivnější?
            #barycenter = np.average(h_plan_sing[Vor_mask_part], axis=0, weights=weights_sim_part)
            h_plan_sing_ma = np.vstack((h_plan_sing_ma, h_plan_sing[Vor_mask_part])) 
            #S_bc = np.cov(h_plan_sing[Vor_mask_part], rowvar=False, bias=True, aweights=weights_sim_part)
            S_bc = np.cov(h_plan_sing_ma, rowvar=False, bias=True, aweights=weights_sim)
            barycenter = np.average(h_plan_sing_ma, axis=0, weights=weights_sim)
            nis_eff += nis
            
        
        
        #print(S_bc)
        #print(nis_ma)
        
        
        
        
        cell_stats = dict()
        # musí sa (na konci) rovnat jedne
        # opravdu dělíme nis'em, jako v normálním IS
        # nikoliv počtem příjatých bodíků, 
        # protože průměrná vaha je o hodně mene významná metrika
        cell_stats['cell_probability'] = np.sum(weights_sim) / nis_eff
        
        
        # tu bacha!
        # takhle se počíta, pokud se netrapíme gradijentem
        # a je to trošiňku optimizovaný, takže čert se nevyzná
        if gradient is None:
            # indexy ii nás moc nezajimajou
            # vzdalenosti snad byjsme zvladli použit?
            
            
            dd2, ii2 = tree.query(h_plan_model_ma, k=[2], p=p_norm)
            dd2 = dd2.flatten()
            ii2 = ii2.flatten()
            
            # tahle hračka s indexy je pro numpy poměrně drahá
            failsii_2 = failsi[ii2]
            # jeden vzorek (včetně hustoty PDF[i])  je nám vždy znám
            # porucha
            if failsi[i]:
                points_1 = PDF[i] * dd2
                node_pf_estimations = points_1 / (points_1 + PDF[ii2] * dd1)
                node_pf_estimations = np.where(failsii_2,1, node_pf_estimations)
                node_pf_pure_estimations = dd2 / (dd1 + dd2)
                node_pf_pure_estimations = np.where(failsii_2,1, node_pf_pure_estimations)
            
                cell_stats['Voronoi_2_point_upper_bound'] = cell_stats['cell_probability'] 
                cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim * node_pf_estimations) / nis_eff
                cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim * node_pf_pure_estimations) / nis_eff
                cell_stats['Voronoi_2_point_lower_bound'] = np.sum(weights_sim[failsii_2]) / nis_eff
                cell_stats['Voronoi_failure_rate'] = cell_stats['cell_probability'] 
                nodes=CandyBox(h_plan.sampling_plan[Vor_mask], w=weights_sim, node_pf_estimations=node_pf_estimations,\
                                    node_pf_pure_estimations=node_pf_pure_estimations, dd1=dd1, dd2=dd2, ii2=ii2)
            
            # neporucha    
            else:
                dd1 = dd1[failsii_2]
                dd2 = dd2[failsii_2]
                points_1 = PDF[i] * dd2
                points_2 = PDF[ii2[failsii_2]] * dd1
                
                node_pf_estimations = points_2 / (points_1 + points_2)
                node_pf_pure_estimations = dd1 / (dd1 + dd2)
            
                cell_stats['Voronoi_2_point_upper_bound'] = np.sum(weights_sim[failsii_2]) / nis_eff
                cell_stats['Voronoi_2_point_failure_rate'] = np.sum(weights_sim[failsii_2]*node_pf_estimations) / nis_eff
                cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(weights_sim[failsii_2] * node_pf_pure_estimations) / nis_eff
                cell_stats['Voronoi_2_point_lower_bound'] = 0
                cell_stats['Voronoi_failure_rate'] = 0
                nodes=CandyBox(h_plan.sampling_plan[Vor_mask][failsii_2], w=weights_sim[failsii_2], node_pf_estimations=node_pf_estimations,\
                                            node_pf_pure_estimations=node_pf_pure_estimations, dd1=dd1, dd2=dd2, ii2=ii2[failsii_2])
                
                # take something with corresponding length
                zeros = np.zeros(len(weights_sim) - len(dd2))
                # add remaining nodes
                nodes.add_sample(CandyBox(h_plan.sampling_plan[Vor_mask][~failsii_2], w=weights_sim[~failsii_2], node_pf_estimations=zeros,\
                                            node_pf_pure_estimations=zeros, ii2=ii2[~failsii_2]))
        
        
        # takhle - pokud chceme gradient použit
        # je třeba eště zoptimalizovať
        else:
            # kolik bodíků jsou nejbližší k mému vzorečkovi
            # np.empty() implicitně má dtype=float
            # tyhle blbosti ponechám jen kvůli callbackovi
            node_pf_estimations = np.empty(len(h_plan_model_ma))
            node_pf_pure_estimations = np.empty(len(h_plan_model_ma))# pure distance estimation
            node_failsi = np.empty(len(h_plan_model_ma), dtype=np.bool) # for L1 Voronoi # co to je za L1 Voronoi?
            
            # projdeme přes každej bodíček
            for node_idx in range(len(h_plan_model_ma)):
                # KDTree byl použit jen k rozdělení na disjunktní úseky, veškerej děj se odehravá tu
                # a to je všechno kvůli gradientu
                node = h_plan_model_ma[node_idx]
                # axis=1 - sčítá všechy směry dohromady, vysledkem je 1D pole rozměru nsim 
                inode2points_model_matrix = np.sum(np.abs(((sampled_plan_model - node) * gradient(node))**p_norm), axis=1)
                #print(inode2points_Rd_matrix)
                
                """
                partition - 
                Creates a copy of the array with its elements rearranged in such a way that
                 the value of the element in k-th position is in the position it would be in a sorted array. 
                 All elements smaller than the k-th element are moved before this element 
                 and all equal or greater are moved behind it. The ordering of the elements in the two partitions is undefined.
                """
                idx = np.argpartition(inode2points_model_matrix, 1) # musí tu bejt 1, počítá sa od nuly
                # je to správný, neboť numpy zaručuje, že druhej prvek (s indexem 1) bude na druhem místě
                node_failsi[node_idx] = failsi[idx[0]] 
                
                                          
                points_weight = PDF[idx[0]] / inode2points_model_matrix[idx[0]] + PDF[idx[1]] / inode2points_model_matrix[idx[1]]
                points_distances = 1 / inode2points_model_matrix[idx[0]] + 1 / inode2points_model_matrix[idx[1]]
                
                failure_weight  = int(failsi[idx[0]]) * PDF[idx[0]] / inode2points_model_matrix[idx[0]]
                failure_weight += int(failsi[idx[1]]) * PDF[idx[1]] / inode2points_model_matrix[idx[1]]
                
                failure_distance = int(failsi[idx[0]]) / inode2points_model_matrix[idx[0]] + int(failsi[idx[1]]) / inode2points_model_matrix[idx[1]]
                
                
                node_pf_estimations[node_idx] = failure_weight/points_weight
                node_pf_pure_estimations[node_idx] = failure_distance/points_distances
            
            
            cell_stats['Voronoi_2_point_upper_bound'] = np.sum(h_plan.w[Vor_mask]*np.ceil(node_pf_estimations)) / nis_eff
            cell_stats['Voronoi_2_point_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_pf_estimations) / nis_eff
            cell_stats['Voronoi_2_point_pure_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_pf_pure_estimations) / nis_eff
            cell_stats['Voronoi_2_point_lower_bound'] = np.sum(h_plan.w[Vor_mask]*np.floor(node_pf_estimations)) / nis_eff
            cell_stats['Voronoi_failure_rate'] = np.sum(h_plan.w[Vor_mask]*node_failsi) / nis_eff
            
            nodes=CandyBox(h_plan.sampling_plan[Vor_mask], w=h_plan.w[Vor_mask], node_pf_estimations=node_pf_estimations,\
                                            node_pf_pure_estimations=node_pf_pure_estimations, node_failsi=node_failsi)
            
            
        
        for key, value in cell_stats.items():
            global_stats[key] += value
        
        # kolbek ↓
        if failsi[i]:
            cell_stats['event'] = 'failure'
        else:
            cell_stats['event'] = 'success'
        callback(estimation=estimation, nodes=nodes, cell_stats=cell_stats, out_nodes=h_plan[~Vor_mask])

        
    
    result = {**estimation, **global_stats}
    # try to put it to the BlackBox
    try:
        sample_box.estimations.append(result)
    except: # asi ani nebudu neuspěch hlasit
        pass
        
    return result
    
    # dedictvi
#for k in range(len(ii)):
#    points_weigths[ii[k]] = points_weigths[ii[k]] + weights_sim[k] / nis
#    near_neighbors[ii[k]] = near_neighbors[ii[k]] + 1
#    Vor_mask[k] = failsi[ii[k]]    
    
    

    
    



def _issi_estimate_outside(f, model_space, sampling_space='G', \
                 powerset_correction=True, design=None, callback=None, outside_budget=1000):
    #č zde f-ko musí taky obsahovat vzorky!
    
    shull = six.Shull(f, model_space, sampling_space, \
                      powerset_correction, incremental=False, design=design)
    
    nodes = shull.integrate(outside_budget)
    
    # -1 = 'out'
    # -2 = 'inside'
    #č převezmeme ISSI
    oiss = shull.oiss
    
    
    # events přepíšeme
    weighted_mean, __, events = oiss.get_means()
    # outside probability
    #č nevyrovnané!
    outside_probability = weighted_mean[events==-1][0]
    
    if callback is not None:
        cell_stats = dict()
        #č tohle je opravdu jen pro formu
        cell_stats['cell_probability'] = outside_probability
        cell_stats['vertex_estimation'] = 0
        cell_stats['weighted_vertex_estimation'] = 0
        cell_stats['event'] = 'outside' #č "outside" se mi libí víc než prostě "out"
        
        
        #:) kolbek ↓
        mask = nodes.is_outside
        #nodes = nodes[mask]
        # out_nodes here meant "not needed"
        out_nodes = nodes[~mask] #č možná je čas se zbavit toho, co se stejně nepouživá?
        callback(shull=shull, nodes=nodes[mask], cell_stats=cell_stats, out_nodes=out_nodes) 
        
    return oiss, outside_probability

    


def full_simplex_estimation(sample_box, model_space='Rn', sampling_space=None, weighting_space=None,\
                           outside_budget=1000, simplex_budget=100, design=None, callback=None):
    """
    Delaunay triangulation
    """
    
    #č zde provadím rozdělení na prostor, ve kterém vzorkujem
    #čs a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
    if sampling_space is None:
        sampling_space = model_space
    
    if weighting_space is None:
        weighting_space = model_space
        
    
    #č jsou to informace (nejen) pro callback 
    estimation={'method': "full_simplex_estimation", 'nsim':sample_box.nsim}
    estimation['model_space'] = model_space
    estimation['sampling_space'] = sampling_space
    estimation['weighting_space'] = weighting_space
    estimation['outside_budget'] = outside_budget
    estimation['simplex_budget'] = simplex_budget
    estimation['design'] = str(design)
    
    
    #č vytahneme ze sample_boxu rozdělení
    f = sample_box.f_model
    
    # powerset_correction=False
    #č totíž my zodpovidáme za to, že Triangulaci pošleme správné ISSI 
    issi, outside = _issi_estimate_outside(f, model_space,\
                     sampling_space=sampling_space, powerset_correction=False,\
                design=design, callback=callback, outside_budget=outside_budget)
                
    tri = six.FullTriangulation(sample_box, tri_space=model_space, issi=issi, \
             weighting_space=weighting_space, incremental=False,\
              on_add_simplex=callback, sampling_space=sampling_space, \
              simplex_budget=simplex_budget, design=design)
              
    tri.integrate()
                
    # TRI-compatible estimation
    # -1 = 'outside', 0=success, 1=failure, 2=mix
    estimations = tri.get_pf_estimation()
    #č doplním
    estimations['global_stats']['outside'] = outside
    
    result = {**estimation, **estimations, 'coplanar':tri.tri.coplanar}
    # try to put it to the BlackBox
    try:
        sample_box.estimations.append(result)
    except: #č asi ani nebudu neuspěch hlasit
        pass
        
    return result
        




def fast_simplex_estimation(sample_box, model_space='Rn', sampling_space=None, weighting_space=None,\
                                            outside_budget=1000, simplex_budget=100, design=None, callback=None):
    
    
    #č vytahneme ze sample_boxu rozdělení
    f = sample_box.f_model
    
    
    #č zde provadím rozdělení na prostor, ve kterém vzorkujem
    #čs a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
    sampled_plan_model = getattr(sample_box, model_space)
    if sampling_space is None:
        sampling_space = model_space
    
    if weighting_space is None:
        weighting_space = model_space
        
    
    #č jsou to informace pro callback 
    estimation={'method': "fast_simplex_estimation", 'nsim':sample_box.nsim}
    estimation['model_space'] = model_space
    estimation['sampling_space'] = sampling_space
    estimation['weighting_space'] = weighting_space
    estimation['outside_budget'] = outside_budget
    estimation['simplex_budget'] = simplex_budget
    estimation['design'] = str(design)
        
    # powerset_correction=True
    #č totíž my zodpovidáme za to, že Triangulaci pošleme správné ISSI 
    issi, outside = _issi_estimate_outside(f, model_space,\
                     sampling_space=sampling_space, powerset_correction=True,\
                design=design, callback=callback, outside_budget=outside_budget)
    
    #č na fastě môžem si dovoliť zkratit cestu 
    failsi = sample_box.failsi
    if np.any(failsi) and not np.all(failsi):
        tri = six.FastTriangulation(sample_box, tri_space=model_space, issi=issi, \
                 weighting_space=weighting_space, incremental=False,\
                  on_add_simplex=callback, sampling_space=sampling_space, \
                  simplex_budget=simplex_budget, design=design)
                  
        tri.integrate()
                    
        # TRI-compatible estimation
        # -1 = 'outside', 0=success, 1=failure, 2=mix
        estimations = tri.get_pf_estimation()
        #č doplním
        estimations['global_stats']['outside'] = outside
        
        result = {**estimation, **estimations, 'coplanar':tri.tri.coplanar}
        
    elif np.all(failsi):
        # TRI-compatible estimation
        # -1=outside, 0=success, 1=failure, 2=mix
        #č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1)
        #č veškerej vnitršek je v poruše
        pf_inside = issi.estimations[-2]
        pf_outside = issi.estimations[-1]
        
        global_stats = dict()
        global_stats['outside'] = outside
        global_stats['success'] = 0
        global_stats['failure'] = 0 # nechce se mi z ISSI vytahovat
        global_stats['mix'] = 0 
        
        # -1=outside, 0=success, 1=failure, 2=mix
        result = {**estimation, \
                'TRI_estimation': {-1:pf_outside, 0:0, 1:pf_inside, 2:0}, \
                'global_stats': global_stats,'vertex_estimation' : pf_inside, \
                'weighted_vertex_estimation' : pf_inside}
        
    else:
        # TRI-compatible estimation
        # -1=outside, 0=success, 1=failure, 2=mix
        #č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1)
        #č vnitršek je asi v pořadku
        pf_inside = issi.estimations[-2]
        pf_outside = issi.estimations[-1]
        
        global_stats = dict()
        global_stats['outside'] = outside
        global_stats['success'] = 0 # nechce se mi z ISSI vytahovat
        global_stats['failure'] = 0
        global_stats['mix'] = 0 
        
        # -1=outside, 0=success, 1=failure, 2=mix
        result = {**estimation, \
                'TRI_estimation': {-1:pf_outside, 0:pf_inside, 1:0, 2:0}, \
                'global_stats': global_stats,'vertex_estimation' : 0, \
                'weighted_vertex_estimation' : 0}
            
    # try to put it to the BlackBox
    try:
        sample_box.estimations.append(result)
    except: #č asi ani nebudu neuspěch hlasit
        pass
        
    return result




#č weighting space nechám pro bodové odhady pf
#č ale podle mně nemá význam
def fast_simplex_cubature(sample_box, tn_scheme, model_space='Rn', sampling_space=None, weighting_space=None,\
                                            outside_budget=1000, callback=None, design=None):
    """
    Delaunay triangulation
    """
    
    #č vytahneme ze sample_boxu rozdělení
    f = sample_box.f_model
    
    
    #č zde provadím rozdělení na prostor, ve kterém vzorkujem
    #čs a prostor "modelu", vô ktôrom, v podstatě, měříme vzdaleností
    sampled_plan_model = getattr(sample_box, model_space)
    if sampling_space is None:
        sampling_space = model_space
    
    if weighting_space is None:
        weighting_space = model_space
        
    
    #č jsou to informace pro callback 
    estimation={'method': "fast_simplex_cubature", 'nsim':sample_box.nsim}
    estimation['tn_scheme'] = tn_scheme.name
    estimation['model_space'] = model_space
    estimation['sampling_space'] = sampling_space
    estimation['weighting_space'] = weighting_space
    estimation['outside_budget'] = outside_budget
    estimation['design'] = str(design)
        
    
    # powerset_correction=True
    #č totíž my zodpovidáme za to, že Triangulaci pošleme správné ISSI 
    issi, outside = _issi_estimate_outside(f, model_space,\
                     sampling_space=sampling_space, powerset_correction=True,\
                design=design, callback=callback, outside_budget=outside_budget)
    
    #č na fastě môžem si dovoliť zkratit cestu 
    failsi = sample_box.failsi
    if np.any(failsi) and not np.all(failsi):
        tri = six.CubatureTriangulation(sample_box, tn_scheme, tri_space=model_space, issi=issi, \
                 weighting_space=weighting_space, incremental=False,\
                  on_add_simplex=callback)
                  
        tri.integrate()
                    
        # TRI-compatible estimation
        # -1 = 'outside', 0=success, 1=failure, 2=mix
        estimations = tri.get_pf_estimation()
        #č doplním
        estimations['global_stats']['outside'] = outside
        
        result = {**estimation, **estimations, 'coplanar':tri.tri.coplanar}
        
    elif np.all(failsi):
        # TRI-compatible estimation
        # -1=outside, 0=success, 1=failure, 2=mix
        #č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1)
        #č veškerej vnitršek je v poruše
        pf_inside = issi.estimations[-2]
        pf_outside = issi.estimations[-1]
        
        global_stats = dict()
        global_stats['outside'] = outside
        global_stats['success'] = 0
        global_stats['failure'] = 0 # nechce se mi z ISSI vytahovat
        global_stats['mix'] = 0 
        
        # -1=outside, 0=success, 1=failure, 2=mix
        result = {**estimation, \
                'TRI_estimation': {-1:pf_outside, 0:0, 1:pf_inside, 2:0}, \
                'global_stats': global_stats,'vertex_estimation' : pf_inside, \
                'weighted_vertex_estimation' : pf_inside}
        
    else:
        # TRI-compatible estimation
        # -1=outside, 0=success, 1=failure, 2=mix
        #č dostaneme vyrovnané odhady Brna-města (-2) a Brna-venkova (-1)
        #č vnitršek je asi v pořadku
        pf_inside = issi.estimations[-2]
        pf_outside = issi.estimations[-1]
        
        global_stats = dict()
        global_stats['outside'] = outside
        global_stats['success'] = 0 # nechce se mi z ISSI vytahovat
        global_stats['failure'] = 0
        global_stats['mix'] = 0 
        
        # -1=outside, 0=success, 1=failure, 2=mix
        result = {**estimation, \
                'TRI_estimation': {-1:pf_outside, 0:pf_inside, 1:0, 2:0}, \
                'global_stats': global_stats,'vertex_estimation' : 0, \
                'weighted_vertex_estimation' : 0}
            
    # try to put it to the BlackBox
    try:
        sample_box.estimations.append(result)
    except: #č asi ani nebudu neuspěch hlasit
        pass
        
    return result