File simplex.py changed (mode: 100644) (index a934b41..8389617) |
... |
... |
from scipy import spatial |
9 |
9 |
from scipy import stats |
from scipy import stats |
10 |
10 |
from .candybox import CandyBox |
from .candybox import CandyBox |
11 |
11 |
|
|
|
12 |
|
import quadpy |
|
13 |
|
|
12 |
14 |
#č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
#č napadlo mě zababáchnuť třidu, která by se sama starala o všem co se tyče |
13 |
15 |
#č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
#č vnější domény. Nešlo mě totíž to udělat jednou funkcí, bylo by velmi |
14 |
16 |
#č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |
#č špatné z hlediska zodpovednosti kódu. Tak to všecko zabalíme to třidy |
|
... |
... |
class _SamplingTriangulation: |
565 |
567 |
|
|
566 |
568 |
|
|
567 |
569 |
class _CubatureTriangulation: |
class _CubatureTriangulation: |
568 |
|
def __init__(sx, sample_box, tn_scheme, tri_space='Rn', issi=None, weighting_space=None, \ |
|
|
570 |
|
def __init__(sx, sample_box, tn_scheme=None, tri_space='Rn', issi=None, weighting_space=None, \ |
569 |
571 |
incremental=True, on_add_simplex=None, on_delete_simplex=None): |
incremental=True, on_add_simplex=None, on_delete_simplex=None): |
570 |
|
|
|
571 |
|
sx.tn_scheme = tn_scheme |
|
572 |
|
|
|
|
572 |
|
|
573 |
573 |
sx.setup(sample_box, tri_space=tri_space, issi=issi,\ |
sx.setup(sample_box, tri_space=tri_space, issi=issi,\ |
574 |
574 |
weighting_space=weighting_space, incremental=incremental,\ |
weighting_space=weighting_space, incremental=incremental,\ |
575 |
575 |
on_add_simplex=on_add_simplex,\ |
on_add_simplex=on_add_simplex,\ |
576 |
576 |
on_delete_simplex=on_delete_simplex) |
on_delete_simplex=on_delete_simplex) |
577 |
|
|
|
578 |
|
|
|
579 |
|
#č vyhodil jsem simplex_id'y |
|
580 |
|
def integrate_simplex(sx, indices, event, event_id, fr, wfr): |
|
581 |
|
#оӵ чылкыт f_model |
|
582 |
|
f = sx.sample_box.f_model |
|
583 |
577 |
|
|
584 |
|
# quadpy |
|
585 |
|
points = [] # container #č vždyť Python nemá pointery |
|
586 |
|
def _get_pdf(x): |
|
587 |
|
fx = f.sample_pdf(x.T, sx.tri_space) |
|
588 |
|
#print("estimate", indices) |
|
589 |
|
#print('x', x.T) |
|
590 |
|
#print('fx', fx) |
|
591 |
|
points.append((x.T, fx)) # side effect |
|
592 |
|
return fx |
|
|
578 |
|
sx.set_tn_scheme(tn_scheme) |
|
579 |
|
|
|
580 |
|
#оӵ чылкыт f_model |
|
581 |
|
#č ptat se na něco u skříňky je extrémně dráho |
|
582 |
|
sx.f = sample_box.f_model |
|
583 |
|
|
|
584 |
|
#č Hodil by se nám levný odhad ve chvili, |
|
585 |
|
#č kdy už si nemôžeme dovoliť víc jak jednoho kandidata v simplexu |
|
586 |
|
#č schema bere vrcholy simplexu, které májí stejné vahy. |
|
587 |
|
#č Tohle by se přece nemohlo selhat! |
|
588 |
|
# tn_fallback_scheme |
|
589 |
|
sx.stroud_tn_1_2 = quadpy.tn.stroud_tn_1_2(sample_box.nvar) |
593 |
590 |
|
|
594 |
|
# premature optimization |
|
595 |
|
if sx.on_add_simplex is not None: |
|
596 |
|
vertices = f[indices] |
|
597 |
|
vertices_model = getattr(vertices, sx.tri_space) |
|
|
591 |
|
|
|
592 |
|
|
|
593 |
|
|
|
594 |
|
|
|
595 |
|
|
|
596 |
|
#č spojení integraci s candidaty nejdřív se zdálo |
|
597 |
|
#č efektivnou opmizací, ale teď začíná obracet protí nám |
|
598 |
|
#č mohli bychom chtit přesně integrovat, |
|
599 |
|
#č ale ukladat jen mírnější počty kandidatů |
|
600 |
|
#č Uděláme to takhle. |
|
601 |
|
def set_tn_scheme(sx, tn_scheme): |
|
602 |
|
sx.tn_scheme = tn_scheme |
|
603 |
|
|
|
604 |
|
if tn_scheme is None: |
|
605 |
|
sx.tn_scheme = sx.stroud_tn_1_2 #č do odhadů se píše tn_scheme.name |
|
606 |
|
sx.integrate_simplex = sx._cheap_integrate_simplex |
|
607 |
|
elif sx.on_add_simplex is not None: |
|
608 |
|
sx.integrate_simplex = sx._callback_integrate_simplex |
598 |
609 |
else: |
else: |
599 |
|
vertices_model = getattr(f, sx.tri_space)[indices] |
|
600 |
|
simplex_measure = sx.tn_scheme.integrate(_get_pdf, vertices_model) |
|
|
610 |
|
sx.integrate_simplex = sx._no_callback_integrate_simplex |
|
611 |
|
|
|
612 |
|
|
|
613 |
|
def _fallback_simplex_integration(sx, vertices): |
|
614 |
|
vertices_model = getattr(vertices, sx.tri_space) |
|
615 |
|
fx = vertices.pdf(sx.tri_space) |
|
616 |
|
|
|
617 |
|
# very special for stroud_tn_1_2 |
|
618 |
|
def _get_pdf(x): return fx |
|
619 |
|
|
|
620 |
|
simplex_measure = sx.stroud_tn_1_2.integrate(_get_pdf, vertices_model) |
601 |
621 |
|
|
602 |
622 |
if (not np.isfinite(simplex_measure)) or (simplex_measure < 0): |
if (not np.isfinite(simplex_measure)) or (simplex_measure < 0): |
603 |
|
print("Kurňa, rozbíla se nám integrace", simplex_measure, points) |
|
|
623 |
|
print("Kurňa, rozbíla se nám integrace totálně", simplex_measure) |
|
624 |
|
print("min_pdf", np.min(fx), "max_pdf", np.max(fx)) |
604 |
625 |
simplex_measure = 0 |
simplex_measure = 0 |
605 |
626 |
|
|
|
627 |
|
return simplex_measure |
|
628 |
|
|
|
629 |
|
# cheap means we will use (existent) simplex vertices |
|
630 |
|
# to integrate by stroud_tn_1_2 cubature scheme |
|
631 |
|
#č skoro není důvod zde použivat quadpy. |
|
632 |
|
#č stroud_tn_1_2 je pouhym průměrem vrcholů |
|
633 |
|
#č ale nechce se mi tahat sem výpočet objemu simplexu s těm determinantem |
|
634 |
|
#č nechame to proto kuadpajovi |
|
635 |
|
def _cheap_integrate_simplex(sx, indices, event, event_id, fr, wfr): |
|
636 |
|
|
|
637 |
|
vertices = sx.f[indices] |
|
638 |
|
simplex_measure = sx._fallback_simplex_integration(vertices) |
|
639 |
|
|
606 |
640 |
assert fr >= 0 |
assert fr >= 0 |
607 |
641 |
assert wfr >= 0 |
assert wfr >= 0 |
608 |
642 |
|
|
|
... |
... |
class _CubatureTriangulation: |
621 |
655 |
cell_stats['event'] = event |
cell_stats['event'] = event |
622 |
656 |
|
|
623 |
657 |
|
|
624 |
|
# kolbek ↓ |
|
625 |
|
nodes_model, pdf = points[0] |
|
626 |
|
nodes = f.new_sample(nodes_model, sx.tri_space) |
|
627 |
|
sx.on_add_simplex(sx, indices=indices, simplex=vertices, nodes=nodes, cell_stats=cell_stats) |
|
|
658 |
|
#č kolbek ↓ |
|
659 |
|
sx.on_add_simplex(sx, indices=indices, simplex=vertices, \ |
|
660 |
|
nodes=vertices, cell_stats=cell_stats) |
|
661 |
|
|
|
662 |
|
|
|
663 |
|
def _callback_integrate_simplex(sx, indices, event, event_id, fr, wfr): |
|
664 |
|
|
|
665 |
|
f = sx.f |
|
666 |
|
|
|
667 |
|
# quadpy |
|
668 |
|
points = [] # container #č vždyť Python nemá pointery |
|
669 |
|
def _get_pdf(x): |
|
670 |
|
nodes = f.new_sample(x.T, sx.tri_space) |
|
671 |
|
fx = nodes.pdf(sx.tri_space) |
|
672 |
|
#print("estimate", indices) |
|
673 |
|
#print('x', x.T) |
|
674 |
|
#print('fx', fx) |
|
675 |
|
points.append((nodes, fx)) # side effect |
|
676 |
|
return fx |
|
677 |
|
|
|
678 |
|
vertices = f[indices] |
|
679 |
|
vertices_model = getattr(vertices, sx.tri_space) |
|
680 |
|
simplex_measure = sx.tn_scheme.integrate(_get_pdf, vertices_model) |
|
681 |
|
|
|
682 |
|
if (not np.isfinite(simplex_measure)) or (simplex_measure < 0): |
|
683 |
|
print("_CubatureTriangulation:") |
|
684 |
|
if simplex_measure < 0: |
|
685 |
|
print("Integráční schema je na houby.", simplex_measure) |
|
686 |
|
else: |
|
687 |
|
print("Kurňa, rozbíla se nám integrace", simplex_measure) |
|
688 |
|
print("min_pdf", np.min(points[0][1]), "max_pdf", np.max(points[0][1])) |
|
689 |
|
simplex_measure = sx._fallback_simplex_integration(vertices) |
|
690 |
|
print("fallback integration result:", simplex_measure) |
|
691 |
|
|
|
692 |
|
assert fr >= 0 |
|
693 |
|
assert wfr >= 0 |
|
694 |
|
|
|
695 |
|
fm = simplex_measure * fr |
|
696 |
|
wfm = simplex_measure * wfr |
|
697 |
|
|
|
698 |
|
#č ISSI tu nemáme, místo toho prostě ukladáme co máme do slovníku |
|
699 |
|
sx.simplex_stats[tuple(indices)] = (event_id, simplex_measure, fm, wfm) |
|
700 |
|
|
|
701 |
|
#č kolbek ↓ |
|
702 |
|
cell_stats = dict() |
|
703 |
|
|
|
704 |
|
cell_stats['cell_probability'] = simplex_measure |
|
705 |
|
cell_stats['vertex_estimation'] = fm |
|
706 |
|
cell_stats['weighted_vertex_estimation'] = wfm |
|
707 |
|
cell_stats['event'] = event |
|
708 |
|
|
|
709 |
|
nodes, pdf = points[0] |
|
710 |
|
sx.on_add_simplex(sx, indices=indices, simplex=vertices, nodes=nodes, cell_stats=cell_stats) |
|
711 |
|
|
|
712 |
|
|
|
713 |
|
|
|
714 |
|
#č spojení integraci s candidaty nejdřív se zdálo |
|
715 |
|
#č efiktivnou opmizací, ale teď začíná obracet protí nám |
|
716 |
|
def _no_callback_integrate_simplex(sx, indices, event, event_id, fr, wfr): |
|
717 |
|
#оӵ чылкыт f_model |
|
718 |
|
f = sx.sample_box.f_model |
|
719 |
|
|
|
720 |
|
# quadpy |
|
721 |
|
def _get_pdf(x): |
|
722 |
|
fx = f.sample_pdf(x.T, sx.tri_space) |
|
723 |
|
return fx |
|
724 |
|
|
|
725 |
|
vertices = f[indices] |
|
726 |
|
vertices_model = getattr(vertices, sx.tri_space) |
|
727 |
|
simplex_measure = sx.tn_scheme.integrate(_get_pdf, vertices_model) |
|
728 |
|
|
|
729 |
|
|
|
730 |
|
if (not np.isfinite(simplex_measure)) or (simplex_measure < 0): |
|
731 |
|
print("_CubatureTriangulation:") |
|
732 |
|
if simplex_measure < 0: |
|
733 |
|
print("Integráční schema je na houby.", simplex_measure) |
|
734 |
|
else: |
|
735 |
|
print("Kurňa, rozbíla se nám integrace", simplex_measure) |
|
736 |
|
simplex_measure = sx._fallback_simplex_integration(vertices) |
|
737 |
|
print("fallback integration result:", simplex_measure) |
|
738 |
|
|
|
739 |
|
|
|
740 |
|
assert fr >= 0 |
|
741 |
|
assert wfr >= 0 |
|
742 |
|
|
|
743 |
|
fm = simplex_measure * fr |
|
744 |
|
wfm = simplex_measure * wfr |
|
745 |
|
|
|
746 |
|
#č ISSI tu nemáme, místo toho prostě ukladáme co máme do slovníku |
|
747 |
|
sx.simplex_stats[tuple(indices)] = (event_id, simplex_measure, fm, wfm) |
628 |
748 |
|
|
629 |
|
|
|
630 |
749 |
|
|
631 |
750 |
|
|
632 |
751 |
|
|