#!/usr/bin/env python
# coding: utf-8

import numpy as np
#from . import IS_stat
from . import sball
#from . import f_models
from scipy import stats
from scipy import spatial
from .candybox import CandyBox


#č jako sůl potřebujem statefull třidu,
#č která by schovávala vnitřní implementaciju,
#č která by nabízela jednotné rozhraní pro ostatní kód WellMet
#č a byla by přiměřeně kompatibilní s ConvexHull ze scipy.spatial
#č Ze scipy bych viděl podporu atributů .points, .npoints, .equations
#č Určitě bych chtěl funkce .update() a .is_outside(sample)
#č Atribute .space by ukazoval, v jakém prostoru konvexní obál je vytvořen


#ё рельса
#č nepodařílo se mi nějak rozumně zobecnit pro libovolný prostor
#č takže Gauss
#č (každá třida implementuje zvlášť)
#def fire(design_point_G, ns):
#    to_fire_G = design_point.G[0]
#    r = np.sqrt(np.sum(np.square(to_fire_G)))
#    a = to_fire_G / r
#    fire_from = stats.norm.cdf(r)
#    t = np.linspace(fire_from, 1, ns)
#    t = stats.norm.ppf(t)
#    fire_G = t.reshape(-1,1) @ a.reshape(1,-1)
#    return design_point.

#č jistě musíme mít nějaký zbytečný kus kódu
#č třida jen pro formu, jen tak na hračku
#č když tečíčky jsou dále jak nejvzdálenější vzorek (bod),
#č tak řekneme, že jsou totálně mimo!
class GBall:
    def __init__(hull, sample):
        hull.sample = sample
        hull.space = 'G' #оӵ малы? Кинлы со кулэ?
    
    @property
    def points(hull):
        return hull.sample.G
        
    @property
    def npoints(hull):
        return len(hull.sample)
        
    @property
    def nsimplex(hull):
        return 1
        
    #def update(hull): pass
    
    def is_inside(hull, nodes):
        R2_hull = np.nanmax(np.sum(np.square(hull.sample.G), axis=1))
        
        R2_nodes = np.sum(np.square(nodes.G), axis=1)
        return R2_nodes < R2_hull 

    def is_outside(hull, nodes): 
        return ~hull.is_inside(nodes)
        
    def get_R(hull):
        return np.sqrt(np.nanmax(np.sum(np.square(hull.sample.G), axis=1)))
        
    def get_r(hull): 
        # zero at best, 
        # on the safe side, between -R and 0
        #r2_hull = np.nanmin(np.sum(np.square(hull.sample.G), axis=1))
        #return -np.sqrt(r2_hull)
        
        #č když na get_r nahlížím z hlediska toho,
        #č že r vymezuje mezikruží, kde nejsme jistí
        #č inside-outside,
        #č tak zde get_r vždy musí vracet prostě R
        return hull.get_R()
        
     
    def fire(hull, ns):
        pass
        
     #č radší nic, než tohle    
#    def fire(hull, ns): # rail
#        R2_hull = np.sum(np.square(hull.sample.G), axis=1)
#        design_index = np.nanargmax(R2_hull)
#        r = np.sqrt(np.nanmax(R2_hull))
#        a = hull.sample.G[design_index] / r
#        fire_from = stats.norm.cdf(r)
#        t = np.linspace(fire_from, 1, ns)
#        t = stats.norm.ppf(t)
#        fire_G = t.reshape(-1,1) @ a.reshape(1,-1)
#        
#        return hull.sample.f_model.new_sample(fire_G, space='G')
        


#č vyhejbám slovu Box, ať nevnáším ještě víc zmatku v označení
class BrickHull: #č nebo BoundingBrick
    def __init__(hull, sample, space='G'):
        hull.sample = sample
        hull.space = space
        
        
        hull._npoints = 0
        hull.mins = np.full(hull.sample.nvar, np.inf)
        hull.maxs = np.full(hull.sample.nvar, -np.inf)
    
    def _update(hull):
        if hull._npoints < hull.npoints:
            hull.mins = np.nanmin(hull.points, axis=0)
            hull.maxs = np.nanmax(hull.points, axis=0)
            hull._npoints = hull.npoints
        
        
    def is_inside(hull, nodes):
        hull._update()
        x = getattr(nodes, hull.space)
        more = np.all(x < hull.maxs, axis=1)
        less = np.all(x > hull.mins, axis=1)
        
        return np.all((more, less), axis=0) #np.all(np.array((more, less)), axis=0)
        
    def is_outside(hull, nodes): 
        return ~hull.is_inside(nodes)
        
    def get_design_points(hull):
        hull._update()
        #sample_model = -hull.A * hull.b.reshape(-1,1)
        sample_model = np.vstack((np.diag(hull.maxs), np.diag(hull.mins)))
        return hull.sample.f_model.new_sample(sample_model, space=hull.space)
        
    # shortcut for Ghull
    # valid only if space==G
    def get_r(hull):
        if hull.space=='G':
            hull._update()
            return np.min((-hull.mins, hull.maxs))
        else:
            return 0
        
    @property
    def points(hull):
        return getattr(hull.sample, hull.space)
        
    @property
    def npoints(hull):
        return len(hull.sample)
        
    @property
    def nsimplex(hull):
        return hull.sample.nvar*2
        
    @property
    def A(hull):
        hull._update()
        #č žádná optimizace, ale kdo tu funkci bude spouštět?
        diag = np.diag(np.ones(hull.sample.nvar))
        return np.vstack((diag, -diag))
        
    @property
    def b(hull):
        hull._update()
        return np.concatenate((hull.maxs, -hull.mins))
        
    @property
    def equations(hull):
        hull._update()
        #č žádná optimizace, ale kdo tu funkci bude spouštět?
        diag = np.diag(np.ones(hull.sample.nvar))
        
        A = np.vstack((diag, -diag))
        b = np.concatenate((-hull.maxs, hull.mins))
        return np.hstack((A,b[:,None]))
        
    def fire(hull, ns): # boom
        sample_U = hull.sample.U
        
        U_mins = np.nanmin(sample_U, axis=0)
        U_maxs = np.nanmax(sample_U, axis=0)
        
        #č min nebo max?
        if np.nanmax(U_mins) > (1-np.nanmin(U_maxs)):
            #č expandujeme se dolu
            var = np.nanargmax(U_mins)
            value = U_mins[var]
            t = np.linspace(value, 0, ns, endpoint=False)
        else:
            #č expandujeme se nahoru
            var = np.nanargmin(U_maxs)
            value = U_maxs[var]
            t = np.linspace(value, 1, ns, endpoint=False)
        
        boom = np.random.random((ns, hull.sample.nvar))
        boom[:, var] = t
        
        return hull.sample.f_model.new_sample(boom, space='U')
    

class DirectHull:
    def __init__(hull, sample, direct_plan, space='G'):
        hull.sample = sample
        hull.direct_plan = direct_plan
        hull.space = space
        
        hull.regen()
        
    def regen(hull):
        hull._npoints = 0
        
        hull.bp = np.full(len(hull.direct_plan), -np.inf)
        #hull.update()
    
    #č nejsem jist, jestli ten update vůbec dělat.
    #č lze navrhnout třidu aj tak, že sama bude hlídat změny.
    #č Jenomže, co kdybychom ten automatický update nechtěli?
    def _update(hull):
        if hull._npoints < hull.npoints:
            #hull.points = getattr(hull.sample, hull.space)
            new_points = hull.points[hull._npoints:]
            
            A = hull.direct_plan @ new_points.T
            new_bp = np.nanmax(A, axis=1)
            hull.bp = np.nanmax((new_bp, hull.bp), axis=0)
            
            hull._npoints = hull.npoints

    @property
    def points(hull):
        return getattr(hull.sample, hull.space)
        
    @property
    def npoints(hull):
        return len(hull.sample)
        
    @property
    def nsimplex(hull):
        return len(hull.direct_plan)
        
    @property
    def A(hull):
        return hull.direct_plan
        
    @property
    def b(hull):
        hull._update()
        return -hull.bp
        
    @property
    def equations(hull):
        return np.hstack((hull.A, hull.b[:,None]))

    def is_inside(hull, nodes):
        hull._update()
        x = getattr(nodes, hull.space)
        
        #E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more)
        A = hull.direct_plan
        bp = np.atleast_2d(hull.bp).T
        
        # N=ns, E - number of hyperplane equations
        ExN = A @ x.T
        higher = np.all(ExN < bp, axis=0)
        return higher
    
    def is_outside(hull, nodes): 
        return ~hull.is_inside(nodes)
        
    def get_design_points(hull):
        hull._update()
        sample_model = -hull.A * hull.b.reshape(-1,1)
        return hull.sample.f_model.new_sample(sample_model, space=hull.space)
        
    # shortcut for Ghull
    # valid only if space==G
    def get_r(hull):
        if hull.space=='G':
            hull._update()
            return np.min(hull.bp)
        else:
            return 0
    
    def fire(hull, ns):
        if hull.space == 'G':
            A = hull.equations[:,:-1]
            b = hull.equations[:,-1]
            
            to_fire = np.nanargmax(b)
            a = A[to_fire]
            fire_from = stats.norm.cdf(hull.get_r())
            t = np.linspace(fire_from, 1, ns, endpoint=False)
            t = stats.norm.ppf(t)
            fire_G = t.reshape(-1,1) @ a.reshape(1,-1)
            
            return hull.sample.f_model.new_sample(fire_G, space='G')


class CompleteHull:
    def __init__(hull, sample, direct_plan, space='G'):
        hull.sample = sample
        hull.direct_plan = direct_plan
        hull.space = space
        
        hull.regen()
        
    def regen(hull):
        hull._npoints = 0
        hull.mins = np.full(hull.sample.nvar, np.inf)
        hull.maxs = np.full(hull.sample.nvar, -np.inf)
        
        hull.bp = np.full(len(hull.direct_plan), -np.inf)
        hull.bm = np.full(len(hull.direct_plan), np.inf)
        
        #hull.update()
    
    #č nejsem jist, jestli ten update vůbec dělat.
    #č lze navrhnout třidu aj tak, že sama bude hlídat změny.
    #č Jenomže, co kdybychom ten automatický update nechtěli?
    def _update(hull):
        if hull._npoints < hull.npoints:
            #hull.points = getattr(hull.sample, hull.space)
            new_points = hull.points[hull._npoints:]
            
            new_mins = np.nanmin(new_points, axis=0)
            new_maxs = np.nanmax(new_points, axis=0)
            
            hull.mins = np.nanmin((new_mins, hull.mins), axis=0)
            hull.maxs = np.nanmax((new_maxs, hull.maxs), axis=0)
            
            A = hull.direct_plan @ new_points.T
            new_bp = np.nanmax(A, axis=1)
            new_bm = np.nanmin(A, axis=1)
            
            hull.bp = np.nanmax((new_bp, hull.bp), axis=0)
            hull.bm = np.nanmin((new_bm, hull.bm), axis=0)
            
            hull._npoints = hull.npoints

    @property
    def points(hull):
        return getattr(hull.sample, hull.space)
        
    @property
    def npoints(hull):
        return len(hull.sample)
        
    @property
    def nsimplex(hull):
        return 2 * (len(hull.direct_plan) + hull.sample.nvar)
        
    @property
    def A(hull):
        hull._update()
        #č žádná optimizace, ale kdo tu funkci bude spouštět?
        diag = np.diag(np.ones(hull.sample.nvar))
        return np.vstack((diag, -diag, -hull.direct_plan, hull.direct_plan))
        
    @property
    def b(hull):
        hull._update()
        return np.concatenate((-hull.maxs, hull.mins, hull.bm, -hull.bp))
        
    @property
    def equations(hull):
        hull._update()
        #č žádná optimizace, ale kdo tu funkci bude spouštět?
        diag = np.diag(np.ones(hull.sample.nvar))
        
        A = np.vstack((diag, -diag, -hull.direct_plan, hull.direct_plan))
        b = np.concatenate((-hull.maxs, hull.mins, hull.bm, -hull.bp))
        return np.hstack((A,b[:,None]))

    def is_inside(hull, nodes):
        hull._update()
        x = getattr(nodes, hull.space)
        
        more = np.all(x < hull.maxs, axis=1)
        less = np.all(x > hull.mins, axis=1)
        
        
        #E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more)
        A = hull.direct_plan
        bp = np.atleast_2d(hull.bp).T
        bm = np.atleast_2d(hull.bm).T
        
        # N=ns, E - number of hyperplane equations
        ExN = A @ x.T
        higher = np.all(ExN < bp, axis=0)
        lower = np.all(ExN > bm, axis=0)
        return np.all((more, less, higher, lower), axis=0)
    
    def is_outside(hull, nodes): 
        return ~hull.is_inside(nodes)
        
    def get_design_points(hull):
        hull._update()
        sample_model = -hull.A * hull.b.reshape(-1,1)
        return hull.sample.f_model.new_sample(sample_model, space=hull.space)
        
    # shortcut for Ghull
    # valid only if space==G
    def get_r(hull):
        if hull.space=='G':
            hull._update()
            return min(-np.max(hull.mins), np.min(hull.maxs), -np.max(hull.bm), np.min(hull.bp))
        else:
            return 0
    
    def fire(hull, ns):
        if hull.space == 'G':
            A = hull.equations[:,:-1]
            b = hull.equations[:,-1]
            
            to_fire = np.nanargmax(b)
            a = A[to_fire]
            fire_from = stats.norm.cdf(hull.get_r())
            t = np.linspace(fire_from, 1, ns, endpoint=False)
            t = stats.norm.ppf(t)
            fire_G = t.reshape(-1,1) @ a.reshape(1,-1)
            
            return hull.sample.f_model.new_sample(fire_G, space='G')
            

class QHull:
    def __init__(self, sample, space='G', incremental=True):
        self.sample = sample
        self.incremental = incremental
        self.space = space
    
    def regen(self):
        points = getattr(self.sample, self.space)
        self.convex_hull = spatial.ConvexHull(points, incremental=self.incremental)
        
    
    def __getattr__(self, attr):
        #č branime se rekurzii
        # defend against recursion
        #оӵ рекурсилы пезьдэт!
        if attr == 'convex_hull':
            #č zkusme rychle konvexní obálky sestavit
            #č a ihned ji vrátit
            self.regen()
            return self.convex_hull
        
        elif attr == 'enough_points':
            return self.sample.nvar < self.sample.nsim
            
        elif attr == 'A':
            return self.convex_hull.equations[:,:-1]
        elif attr == 'b':
            return self.convex_hull.equations[:,-1]
            
        elif attr == 'points':
            if self.enough_points:
                return self.convex_hull.points
            else:
                return getattr(self.sample, self.space)
                
        elif attr == 'npoints':
            if self.enough_points:
                return self.convex_hull.npoints
            else:
                return len(self.sample)
        
        elif attr == 'nsimplex':
            if self.enough_points:
                return self.convex_hull.nsimplex
            else:
                return 0
            
            
        #ё По всем вопросам обращайтесь 
        #ё на нашу горячую линию    
        else:
            self._update() #č dycky čerstý chleba!
            return getattr(self.convex_hull, attr)
            
    
    #č nejsem jist, jestli ten update vůbec dělat.
    #č lze navrhnout třidu aj tak, že sama bude hlídat změny.
    #č Jenomže, co kdybychom ten automatický update nechtěli?
    def _update(self):
        if self.convex_hull.npoints < self.sample.nsim:
            if self.incremental:
                points = getattr(self.sample, self.space)  
                self.convex_hull.add_points(points[self.convex_hull.npoints:])
            else:
                self.regen()
        
        
    def is_inside(self, nodes):
        if self.enough_points:
            self._update()
            x = getattr(nodes, self.space)
            
            #E [normal, offset] forming the hyperplane equation of the facet (see Qhull documentation for more)
            A = self.convex_hull.equations[:,:-1]
            b = self.convex_hull.equations[:,-1]
            
            # N=ns, E - number of hyperplane equations
            ExN = A @ x.T + np.atleast_2d(b).T
            mask = np.all(ExN < 0, axis=0)
            return mask
        else:
            return np.full(len(nodes), False)
    
    def is_outside(hull, nodes): 
        return ~hull.is_inside(nodes)
        
    def get_design_points(hull):
        hull._update()
        sample_model = -hull.A * hull.b.reshape(-1,1)
        return hull.sample.f_model.new_sample(sample_model, space=hull.space)
        
    def get_design_point(hull):
        hull._update()
        to_fire = np.nanargmax(hull.b) #č tak, pro jistotu
        sample_model = -hull.A[to_fire] * hull.b[to_fire]
        return hull.sample.f_model.new_sample(sample_model, space=hull.space)
    
        
    # shortcut for Ghull
    # valid only if space==G
    def get_r(hull):
        if hull.space=='G':
            if hull.enough_points:
                hull._update()
                b = hull.convex_hull.equations[:,-1]
                return -np.nanmax(b)
            else:
                return 0
        else:
            return 0
    
    def fire(hull, ns):
        if hull.space == 'G':
            try:
                A = hull.equations[:,:-1]
                b = hull.equations[:,-1]
                
                to_fire = np.nanargmax(b) #č tak, pro jistotu
                a = A[to_fire]
                fire_from = stats.norm.cdf(hull.get_r())
                t = np.linspace(fire_from, 1, ns, endpoint=False)
                t = stats.norm.ppf(t)
                fire_G = t.reshape(-1,1) @ a.reshape(1,-1)
                
                return hull.sample.f_model.new_sample(fire_G, space='G')
            except:
                pass


#č mým úkolem při návrhu této třidy je pořádně všecko zkomplikovat.
#č Dostávám za to peníze. 
#č Takže. Udělám 3 druhů estimátorů
# convex_hull_estimation  -2: inside, -1: outside
# shell_estimation  -22: inner, -3: shell, -11: outer
# ghull_estimation  -22: inner, -21: shell inside, -12: shell outside, -11: outer 

try: # try to outthink f****** OS's memory management
    import os
    mem_bytes = os.sysconf('SC_PAGE_SIZE') * os.sysconf('SC_PHYS_PAGES')
except BaseException as e:
    mem_GB = 16
    mem_bytes = mem_GB * 1024**3 # hello, Windows!
    print("convex_hull failed to get amount of RAM installed. %s GB is supposed."% mem_GB, repr(e))
    print("BTW, you are using Windows, aren't you?")
    

class Ghull:
    def __init__(self, hull, non_Gaussian_reduction=0.9, design=None):
        self.hull = hull
        self.design = design
        self.shell = sball.Shell(hull.sample.nvar)
        self.outside_dist = sball.Shell(hull.sample.nvar)
        self.sample = hull.sample
        
        self._r = 0 # estimation for nonGaussian distributions
        self.non_Gaussian_reduction = non_Gaussian_reduction
        self.integration_cutoff = np.inf
        
    def fire(self, ns):
        nodes = self.hull.fire(ns)
        if nodes is not None:
            return nodes
        else:
            return self.boom(ns)
            
    def boom(self, ns):
        self.outside_dist.set_bounds(self.get_R())
        nodes_G = self.outside_dist.rvs(ns)
        nodes = self.sample.f_model.new_sample(nodes_G, space='G')
        return nodes
        
    def get_R(self):
        sum_squared = np.sum(np.square(self.sample.G), axis=1)
        #index = np.argmax(sum_squared)
        return np.sqrt(np.nanmax(sum_squared))
        
    def get_r(self):
        "calculates minimal signed distance to planes. Can therefore be negative"
        #return -np.nanmax(self.hull.b)
        if self.hull.space == 'G':
            return self.hull.get_r()
        else:
            return self._r * self.non_Gaussian_reduction
     
    def get_shell_estimation(self):
        shell = self.shell
        r = self.get_r()
        R = self.get_R()
        
        if r<0:
            shell.set_bounds(0, R)
        else:
            shell.set_bounds(r, R)
        
        # shell_estimation  -22: inner, -3: shell, -11: outer
        shell_estimation = {-22:shell.ps, -3: shell.p_shell, -11: shell.pf}
        global_stats = {"nsim":self.sample.nsim, "ndim":self.sample.nvar, \
                        "nfacets": self.hull.nsimplex, "r":r, "R":R, \
                    "inner":shell.ps, "shell":shell.p_shell, "outer":shell.pf}
        return shell_estimation, global_stats
        
    def integrate(self, nis):
        candy_nodes, _nis, nf, shell_estimation, global_stats = self.check_it_out(nis)
        assert _nis == nis
        
        ns = nis - nf
        # shell_estimation  -22: inner, -3: shell, -11: outer
        p_shell = shell_estimation.pop(-3)
        shell_pf = nf/nis * p_shell
        shell_ps = ns/nis * p_shell
        
        # ghull_estimation  -22: inner, -21: shell inside, -12: shell outside, -11: outer 
        ghull_estimation = shell_estimation; del(shell_estimation)
        ghull_estimation[-21] = shell_ps
        ghull_estimation[-12] = shell_pf
        global_stats["shell_budget"] = nis
        global_stats["shell_inside"] = shell_ps
        global_stats["shell_outside"] = shell_pf
        
        # convex_hull_estimation  -2: inside, -1: outside
        inside = shell_ps + self.shell.ps
        outside = shell_pf + self.shell.pf
        convex_hull_estimation = {-2: inside, -1: outside}
        global_stats["inside"] = inside
        global_stats["outside"] = outside
        
        return candy_nodes, ghull_estimation, convex_hull_estimation, global_stats
        
        
    def rvs(self, nsampled, seats, ns): 
        "Generování­ vzorků (kandidátů a integračních bodů)"
        # rand_dir: prepare ns random directions on a unit d-sphere
        rand_dir = sball.get_random_directions(seats, self.sample.nvar) #random directions
        
        # generate sampling probabilites
        left = (ns-nsampled) / ns
        right = (ns-nsampled-seats) / ns
        #č přidáme trochu zmatku.
        #č globálně jdeme r->R, localně ale R_i+ -> R_i-
        #č menší poloměry musejí jít dřív - na to zavazano nonGaussian _r
        #č převracení lokalně umožní vždycky mít alespoň jeden bodík outside,
        #č což je taky velmi vhodné, vzhledem k tomu, že _r bere z outside bodíků
        p = np.linspace(right, left, seats, endpoint=False) # probabilities for the radius
        
        # convert probabilitites into random radii
        # (distances from origin that are greater than r and less than R)
        r = self.shell.isf(p) # actually, it is the same as CDF inverse
        
        #finally a random sample from the optimal IS density:
        sample_G = rand_dir*r[:,None]
        
        return sample_G
        
    def check_it_out(self, nis):
        #č no to teda disajn třidy je doopravdy hroznej
        # .get_shell_estimation() will calculate radiuses and will update shell
        shell_estimation, global_stats = self.get_shell_estimation()
        
        #č rvs má vzorkovat od měnšího poloměru k většímu.
        #č to znamená, že když první outside bude mít nejměnší r vůbec
        r_needed = (self.hull.space!='G')
        
        
        nsampled = 0
        nfailed = 0
        if self.hull.nsimplex == 0:
            bus = self.integration_cutoff
        else:
            bus = int(mem_bytes / self.hull.nsimplex / 8 / 3) + 1
        while nsampled < nis:
            
            seats = min(nis - nsampled, self.integration_cutoff, bus)
            try:
                nodes_G = self.rvs(nsampled, seats, nis)
                nodes = self.sample.f_model.new_sample(nodes_G, space='G')
                mask = self.hull.is_outside(nodes)
                if r_needed and np.any(mask):
                    self._r = np.sqrt(np.nanmax(np.sum(np.square(nodes_G[mask]), axis=1)))
                    r_needed = False
                # -2 = 'inside' -1 = 'outside'
                candy_nodes = CandyBox(nodes, event_id=mask-2, is_outside=mask)
                
                #č nevím proč, ale kdysi mě to vyšlo rychlejc jak obyčejný součet
                nfailed += len(mask[mask]) 
                nsampled += len(mask) 
                assert len(mask) == seats
            except MemoryError as m:
                print("Ghull: memory error, %s sedaček" % seats, repr(m))
                self.integration_cutoff = int(np.ceil(seats/3))
        
        #č Candy bodů nemusí být nis!
        return candy_nodes, nsampled, nfailed, shell_estimation, global_stats


##č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče 
##č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi
##č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy  
##č a odebereme z už beztak přetíženého blackboxu část komplexity
## keywords: ISSI, estimation, outside, ConvexHull, Sball, IS kolem středních hodnot
#class Shull: # issi_estimate_outside
#    def __init__(sx, f, model_space, sampling_space='G', \
#                 powerset_correction=True, incremental=True, design=None):
#        #č tím powerset_corretion je myšlena úplná soustava jevů,
#        #č tj. vyrovnaní s použitím míry vnější i vnitřní
#        #č powerset_correction=True přídá -2 (inside) jev do ISSI
#        
#        #č zde f-ko musí taky obsahovat vzorky!
#        sx.f = f
#        sx.model_space = model_space
#        
#        #č kašlu na to, pokud uživatel zadal nesmysl, tak
#        #č sam převedu do nečěho smyslúplnějšího
#        _dict = {'R':'Rn', 'aR':'aRn', 'P':'GK', 'aP':'aGK', 'U':'G', 'aU':'aG'}
#        if sampling_space in _dict:
#            sx.sampling_space = _dict[sampling_space]
#        else:
#            sx.sampling_space = sampling_space
#            
#        sx.design = design
#            
#        sampled_plan_model = getattr(f, model_space)
#        #č žádná kontrola chyb - nechť to spadné, když má spadnout!
#        sx.convex_hull = spatial.ConvexHull(sampled_plan_model, incremental=incremental)
#        
#        # current outside probability estimation
#        sx.p_out = 0.5 # implicit value
#        sx.sball = sball.Sball(f.nvar)
#        sx.base_r = sx.sball.get_r(0.5) # we want in average 50/50 ratio
#        
#            
#        # -1 = 'outside'
#        # -2 = 'inside'
#        #č založme ISSI
#        sx.powerset_correction = powerset_correction
#        #č potřebuji pro korektnost mít před integrací zadané jevy 
#        if powerset_correction:
#            sx.oiss = IS_stat.ISSI([-1, -2]) 
#        else:
#            sx.oiss = IS_stat.ISSI([-1]) 
#        
#        
#    
#    def increment(sx, input_sample):
#        #č sample by měl byt jíž převeden na f (v .add_sample()),
#        #č zodpovidá za to volajicí kód!
#        sx.convex_hull.add_points(getattr(input_sample, sx.model_space))
#        
#        
#        
#        # require
#    def integrate(sx, nis):
#        # getting rid of the old estimations
#        sx.oiss.delete_event_data(-1)
#        sx.oiss.events.append(-1) #č už tak trošku sahám do vnitřku cizí třidy
#        if sx.powerset_correction:
#            sx.oiss.delete_event_data(-2)
#            sx.oiss.events.append(-2) #č a záse
#    
#        #č posunutí středu trošiňku porušuje předpoklady, za kterých 
#        #č Sball volí rozdělení, ale přečo se mi stavá,
#        #č že ve více dimenzích Shull na začatku prostě nemůže 
#        #č trefit ten blbej... v podstatě simplex
#        #č Těžiště, přesnějí, prostě nějaký střed můžeme najít dvěma způsoby: 
#        #č 1. mean(vertices.model_space).sampling_space
#        #č 2. mean(vertices.sampling_space)
#        #č myslím si, že ten první (a bez váh) je stabilnější 
#        #č (víme, jak střed vrcholů v nějakém prostoru může ani netrefit do 
#        #č simplexu v původním prostoru)
##        vertices_model = sx.convex_hull.points[sx.convex_hull.vertices]
##        barycenter_model = np.mean(vertices_model, axis=0)
##        if sx.model_space == sx.sampling_space:
##            sx.barycenter_sampling = barycenter_model
##        else:
##            barycenter = sx.f.new_sample(barycenter_model, sx.model_space)
##            sx.barycenter_sampling = np.squeeze(getattr(barycenter, sx.sampling_space))
#
#        #č rozhodl jsem, že shull musí odhadovat outside, ne motat se centrem Brna.
#        #č předpokladám, že uživatel zadal buď G, nebo aspoň Rn prostor
#        #sx.barycenter_sampling = np.full(sx.f.nvar, 0, dtype=np.int8)
#        
#        # first step
#        nodes = sx._sample_sball(nis)
#        mask = nodes.is_outside
#        
#        
#        cut_off_out = int(nis/3)
#        cut_off_in = int(nis/3)
#        #č robím cyklus dokud nesberu dostatečně teček.
#        #č To je fakt nejrobustnější řešení, co mě napadá
#        # while (number_of_out_nodes or number_of_nodes_inside is too_small)
#        while (len(mask[mask]) < cut_off_out): # or (len(mask[~mask]) < cut_off_in):
#            print(sx.__class__.__name__ + ":", "cut off! Outsides: %s, insides: %s, p_out=%s"\
#                                                 % (len(mask[mask]), len(mask[~mask]), sx.p_out))
#            #č je třeba jenom sehnat dostatečně bodíků a utikat
#            nodes.add_sample(sx._sample_sball(nis))
#            mask = nodes.is_outside
#            
#        #č když neprovadíme výrovnání, tak ten vnitršek nachren nepotřebujem
#        if (len(mask[~mask]) < cut_off_in) and sx.powerset_correction:
#            print(sx.__class__.__name__ + ":", \
#            "cut off inside (%s of %s needed)! Do simplex-based integration of convex hull"\
#                                                     % (len(mask[~mask]), cut_off_in))
#            
#            nodes.add_sample(sx._sample_simplex(nis))
#        
#            
#        return nodes
#        
#        
#        
#    def _sample_simplex(sx, nis):
#        
#        #č f-ko sice musí odkazovat na aktuální f_model
#        #č ale na druhou stranu normálně ._integrate_simplex()
#        #č potřebujeme pouze jednou - hned při vytvaření
#        vertices = sx.f[sx.convex_hull.vertices]
#        nvar = vertices.nvar
#        
#        # IS_like uses .new_sample method, so vertices can not be a SampleBox object
#        #
#        #č IS_like těžiště počítá v sampling_space, ale asi mu to až tak nevadí
#        #
#        # already divided by nsim in variance formule
#        # divide by /(nvar+1)/(nvar+2) from simplex inertia tensor solution
#        # multiply by simplex_volume, but it looks like it shouldn't be here
#        # for simplex: d = nvar+2 
#        #č sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku
#        nodes = IS_stat.IS_like(vertices, sampling_space=sx.sampling_space, \
#                nis=nis, d=nvar+2, design=sx.design)
#        
#        #č indikatorová funkce
#        sx.is_outside(nodes)
#        
#        # for IS_stats
#        #č zkoušel jsem zadavat celou serii - zhoršuje to odhady
#        #č nemůžeme důvěrovat tomu, jak ten slepej simplex vidí svět
#        weights = nodes.w[~nodes.is_outside]
#        sx.oiss.add_single_event_data(weights, event=-2, nis=nis)
#        # add_single_event_data() do not calculate estimations itself
#        sx.oiss.get_estimations()
#        return nodes
#    
#    
#    
#    def _sample_sball(sx, nis):
#        nvar = sx.f.nvar
#        sampling_r, sx.p_out = sx.sball.get_r_iteration(sx.p_out)
#        #sampling_r = sx.sball.get_r(sx.p_out)
#        
#        #č asi tam bylo sampling_r/base_r, že?
#        std_ball = sampling_r/sx.base_r
#        #č chcu std=1, když p_out -> 1
#        #č a std=sball_solušn, když p_out -> 0
#        #č surovější nevymyslíš! 
#        std = std_ball + sx.p_out
#        #č u stats.norm zadáváme směrodatnou odchylku, je to asi správné
#        #h = f_models.UnCorD([stats.norm(sx.barycenter_sampling[i], std) for i in range(nvar)])
#        #nodes = IS_stat.IS(sx.f, h, space_from_h='R', space_to_f=sx.sampling_space,  Nsim=nis)
#        
#        #norm_params = [(sx.barycenter_sampling[i], std) for i in range(nvar)]
#        nodes = IS_stat.IS_norm(sx.f, std=std, \
#                sampling_space=sx.sampling_space, nis=nis, design=sx.design)
#        
#        outside_measure = sx._apply_nodes(nodes, nis)
#        
#        #č pro přiště
#        sx.p_out = (sx.p_out + outside_measure) / 2
#        
#        return nodes
#    
#    
#    def _apply_nodes(sx, nodes, nis):
#        #č indikatorová funkce
#        sx.is_outside(nodes)
#        
#        # for IS_stats
#        if sx.powerset_correction:
#            #č získáme výrovnaný odhad - je to skoro zdarma
#            #svar = (sampling_r/sx.base_r)**2 # svar like sampling_variance
#            #č kdysi snažil jsem něco odvést, moc se mi to nepovedlo
#            #č je to jen tak, jeden z pokusu, hrubej nastřel
#            #im = svar**nvar * np.exp(nvar/svar - nvar)
#            
#            #č radší ne. IM špatně zachycuje nizkou důvěru k tomu, co nemá vlastní tečky
#            sx.oiss.add_IS_serie(nodes.w, nodes.event_id, implicit_multiplicator=np.inf)
#            outside_measure = sx.oiss.estimations[-1]
#        else:
#            weights = nodes.w[nodes.is_outside]
#            #č IM všecko pokazí, jakmile začnu přídávát další jevy
#            sx.oiss.add_single_event_data(weights, event=-1, nis=nis)
#            # add_single_event_data() do not calculate estimations itself
#            weighted_mean, __, events = sx.oiss.get_means()
#            # outside probability
#            #č nevyrovnané!
#            outside_measure = weighted_mean[events==-1][0]
#            
#        return outside_measure
#    
#    
#    
#    def is_outside(sx, nodes):
#        node_coordinates = getattr(nodes, sx.model_space)
#        mask = is_outside(sx.convex_hull, node_coordinates)
#        
#        # -1 = 'outside'
#        # -2 = 'inside'
#        nodes.event_id = mask - 2
#        nodes.is_outside = mask
#
#