iam-git / WellMet (public) (License: MIT) (since 2021-08-31) (hash sha1)
WellMet is pure Python framework for spatial structural reliability analysis. Or, more specifically, for "failure probability estimation and detection of failure surfaces by adaptive sequential decomposition of the design domain".
Clone URLs: https://rocketgit.com/user/iam-git/WellMet ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet git://git.rocketgit.com/user/iam-git/WellMet
release trunk
List of commits:
Subject Hash Author Date (UTC)
convex_hull.QHull: implement erfc and gammaincc-based outside integrations for now ad0a665b53b6c8ca24a215c0e4ca73a5c1d145ad I am 2023-01-23 19:29:09
convex_hull.QHull: accept integrational crime 788f3ffe604a4cac9d19382250dd5ebe3cd5befc I am 2023-01-21 13:09:31
convex_hull.QHull: normal fix 083bbdc4cb41fdb3e357b1d0359705a2a65d6d27 I am 2023-01-20 19:07:17
convex_hull.QHull: implement chi cubature integration 7537b642967b7dd5d8475c20aff2a930a5bcd3ad I am 2023-01-20 15:54:05
convex_hull: implement QHullCubature class 6f1b51ba6fe281f7cae4f44889db96787d916a99 I am 2023-01-20 10:46:54
convex_hull.QHull: add proof-of-concept get_inside() method 6013fbac742d6070771d7798c928538ebc54495e I am 2023-01-19 10:19:32
whitebox: add r exact for common nD boxes 7c4ab62f069f0efeb48a403dac5d27033003ded2 I am 2023-01-17 15:09:08
qt_gui: improve GRaph to show r-R bounds 59fa5bde7be2197b1f739bd1f829f927cef62603 I am 2023-01-17 14:26:16
dicebox.circumtri: calculate event's radia c6868cfbde7c8a86b73217da544f60e94a117bd8 I am 2023-01-16 06:12:34
simplex.Triangulation.get_averaged_mixed_normals: implement deep averaging 6dfb056e477479331c7a903ea52a98301fca4e71 I am 2023-01-15 23:26:43
replace flatten() by reshape(-1) everywhere. The former always allocates new array 97b35c989a5ad8e459bb60133222d50c7f78e1f0 I am 2023-01-15 12:04:30
qt_gui.qt_gui: explicitly export sample_box to console's namespace dd2e50245157082ddb645b693a8a8c573215d50c I am 2023-01-15 11:42:01
replace everywhere np.bool by just bool (it looks like newer numpy dropped it) ac208f89bbf0ad76d38790686041ea778a9b6376 I am 2023-01-15 11:31:56
simplex._Triangulation: implement averaged gradients() f6608597cf490dd6cf1a6af3b7d8e6d0a72732d9 I am 2023-01-15 09:37:14
add line and two lines problems for sensitivity analysis 0c620f7960399ed8b34e02d663772d0eb34e66f1 I am 2023-01-14 08:30:56
simplex.FullCubatureIntegration.get_sensitivities: do not normalize global gradient c648c484a5599a00236ab70c36cdf407c24f5306 I am 2023-01-14 04:51:16
simplex: implement sensitivities 65056319c51fd93deb76ce14e52c08869dfc04a2 I am 2023-01-13 11:27:46
mplot.maxes: prepare GRaph plot 68774535915af3deaf6104a0bbe9b3c4d24c0cec I am 2023-01-12 02:43:20
mplot.mgraph.tri_estimation_plot: use planar vertex estimator 1d300fa56248be9cdab5a4e5b0d4372f68fd5255 I am 2023-01-11 15:48:28
mplot.mart: implement lsf_boundary() contouring 8642d2f0c043f4e71af1cc45b9e25da326afd3ca I am 2023-01-11 05:20:00
Commit ad0a665b53b6c8ca24a215c0e4ca73a5c1d145ad - convex_hull.QHull: implement erfc and gammaincc-based outside integrations for now
Author: I am
Author date (UTC): 2023-01-23 19:29
Committer name: I am
Committer date (UTC): 2023-01-23 19:29
Parent(s): 788f3ffe604a4cac9d19382250dd5ebe3cd5befc
Signer:
Signing key:
Signing status: N
Tree: 19b5be46373aab8a98507bd2a6927a6828d28216
File Lines added Lines deleted
wellmet/convex_hull.py 111 33
File wellmet/convex_hull.py changed (mode: 100644) (index d23aee9..d12f786)
... ... import numpy as np
5 5 import mpmath import mpmath
6 6 from scipy import stats from scipy import stats
7 7 from scipy import spatial # for QHull from scipy import spatial # for QHull
8 from scipy import special # for QHull outside integration
8 9
9 10 from . import sball from . import sball
10 11
 
... ... QHullEstimation = namedtuple('QHullEstimation', (
825 826 )) ))
826 827
827 828
829 def get_nsphere_const(nvar, const=1):
830 if nvar == 2:
831 return const * 2 * np.pi
832 elif nvar == 1:
833 return const * 2
834 else:
835 return get_nsphere_const(nvar - 2, const * 2 * np.pi / (nvar - 2) )
836
828 837 class QHull: class QHull:
829 838 # take some global function # take some global function
830 839 _orth_helper = _orth_helper _orth_helper = _orth_helper
 
... ... class QHull:
840 849 self.incremental = incremental self.incremental = incremental
841 850 self.space = space self.space = space
842 851 self.auto_update = auto_update self.auto_update = auto_update
843 self.tn_scheme = tn_scheme
852
844 853 self.sball = sball.Sball(sample.nvar) self.sball = sball.Sball(sample.nvar)
845 854 self.shell = sball.Shell(sample.nvar) self.shell = sball.Shell(sample.nvar)
846 855 self.fallback_plan = quadpy.un.mysovskikh_1(sample.nvar).points self.fallback_plan = quadpy.un.mysovskikh_1(sample.nvar).points
856
857 if tn_scheme is None:
858 self.tn_scheme = quadpy.tn.grundmann_moeller(sample.nvar - 1, 5)
859 else:
860 self.tn_scheme = tn_scheme
861
862 nvar = sample.nvar
863 #self.chi_const = 2**((nvar - 1) / 2) * special.gamma(nvar - 1/2) / special.gamma(nvar / 2)
864 self.chi_like_const = np.sqrt(np.pi) / np.sqrt(2) / (2 * np.pi)**(nvar/2)
865 self.nsphere_surface_area = get_nsphere_const(nvar)
866
867
847 868
848 869 def regen(self): def regen(self):
849 870 points = getattr(self.sample, self.space) points = getattr(self.sample, self.space)
 
... ... class QHull:
1028 1049 #č neřikejte mi, že musím pěčlivějc zpracovavat chyby! #č neřikejte mi, že musím pěčlivějc zpracovavat chyby!
1029 1050 return 1 return 1
1030 1051
1031 def get_norm_inside(hull):
1052 def get_norm_outside(hull):
1032 1053 hull._update() hull._update()
1033 1054 if hull.space != 'G': if hull.space != 'G':
1034 1055 raise raise
1035 1056
1036 1057 vertices_model = hull.points[hull.simplices].transpose((1, 0, 2)) vertices_model = hull.points[hull.simplices].transpose((1, 0, 2))
1037 result = hull.tn_scheme.integrate(hull._norm_inside_callback, vertices_model)
1058 result = hull.tn_scheme.integrate(hull._outside_callback, vertices_model)
1038 1059
1039 return np.sum(result, axis=1)
1060 return np.sum(result) * hull.chi_like_const
1040 1061
1041 1062
1042 1063 def get_chi_outside(hull): def get_chi_outside(hull):
1043 1064 hull._update() hull._update()
1044 1065 if hull.space != 'G': if hull.space != 'G':
1045 1066 raise raise
1046 scheme = quadpy.tn.grundmann_moeller(hull.sample.nvar - 1, 5)
1067
1047 1068
1048 1069 vertices_model = hull.points[hull.simplices].transpose((1, 0, 2)) vertices_model = hull.points[hull.simplices].transpose((1, 0, 2))
1049 1070 result = hull.tn_scheme.integrate(hull._chi_outside_callback, vertices_model) result = hull.tn_scheme.integrate(hull._chi_outside_callback, vertices_model)
1050 1071
1051 return np.sum(result)
1072 return np.sum(result) / hull.nsphere_surface_area
1073
1074 # # quadpy
1075 # def _norm_inside_callback(hull, x):
1076 # # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56)
1077 #
1078 # cdfs = stats.norm.cdf(x)
1079 # pdfs = stats.norm.pdf(x)
1080 #
1081 # for i in range(hull.sample.nvar):
1082 # cdfs[i] *= np.prod(pdfs[:i], axis=0)
1083 # cdfs[i] *= np.prod(pdfs[i+1:], axis=0)
1084 # cdfs[i] *= np.atleast_2d(hull.A[:, i]).T
1085 #
1086 #
1087 # # n_values x nsimplex x scheme_points
1088 # return cdfs
1089
1052 1090
1053 # quadpy
1054 def _norm_inside_callback(hull, x):
1091
1092
1093 def _outside_callback(hull, x):
1055 1094 # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56) # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56)
1056 1095
1057 cdfs = stats.norm.cdf(x)
1058 pdfs = stats.norm.pdf(x)
1096 #č přijde velké, obrovské numpy pole
1097 #č a já bych chtěl, aby se tu zbytečně nealokovalo
1098
1099 result = np.einsum('ijk,ijk->jk', x, x) # r**2
1100 np.sqrt(result, out=result) # r
1059 1101
1060 for i in range(hull.sample.nvar):
1061 cdfs[i] *= np.prod(pdfs[:i], axis=0)
1062 cdfs[i] *= np.prod(pdfs[i+1:], axis=0)
1063 cdfs[i] *= np.atleast_2d(hull.A[:, i]).T
1102 #result = r / np.sqrt(2) # r / sqrt(2)
1103 np.divide(result, np.sqrt(2), out=result)
1104 #special.gammaincc(hull.sample.nvar - 1/2, result, out=result) # F
1105 special.erfc(result, out=result)
1064 1106
1107 # #č normalizujeme normály. Jenže normály samotné dodáme později
1108 # result /= r
1109 #
1110 # # A = nsimplex x nvar
1111 # # x = nvar x nsimplex x scheme_npoints
1112 # # normals = nsimplex x scheme_npoints
1113 # #normals = np.einsum('ij,jik->ik', hull.A, x)
1114 # np.einsum('ij,jik->ik', hull.A, x, out=r) # normals
1115 # result *= r
1065 1116
1066 1117 # n_values x nsimplex x scheme_points # n_values x nsimplex x scheme_points
1067 return cdfs
1118 return result
1068 1119
1069
1070 # quadpy
1071 def _chi_outside_callback(hull, x):
1120 def _outside_normals_callback(hull, x):
1072 1121 # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56) # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56)
1073 1122
1074 #č zde je zločin:
1075 #č r-ko by se mělo dostat do integralu k chi_pdf
1076 #č a místo survival function chi rozdělení by se mělo
1077 #č nějaká jiná primitivní funkce
1078 r = np.sqrt(np.sum(np.square(x), axis=0))
1079 cdfs = hull.sball.sf(r)
1080 # chi_pdf * bůhvíco == phi(x) * phi(y) * phi(z)
1081 # bůhvíco == 2 * pi * r
1082 #chi_pdf = stats.chi.pdf(r, nvar)
1083 #norm_pdf = np.prod(stats.norm.pdf(x), axis=0)
1123 #č přijde velké, obrovské numpy pole
1124 #č a já bych chtěl, aby se tu zbytečně nealokovalo
1084 1125
1085 r_vectors = x / r[None, :, :]
1126 r = np.einsum('ijk,ijk->jk', x, x) # r**2
1127 np.sqrt(r, out=r) # r
1086 1128
1087 normals = np.sum(r_vectors * hull.A.T[:, :, None], axis=(0))
1129 result = r / np.sqrt(2) # r / sqrt(2)
1130 #np.divide(result, np.sqrt(2), out=result)
1131 #special.gammaincc(hull.sample.nvar - 1/2, result, out=result) # F
1132 special.erfc(result, out=result)
1088 1133
1089 # print(norm_pdf / chi_pdf)
1090 # print(0.5 / r / np.pi)
1134 #č normalizujeme normály. Jenže normály samotné dodáme později
1135 result /= r
1136
1137 # A = nsimplex x nvar
1138 # x = nvar x nsimplex x scheme_npoints
1139 # normals = nsimplex x scheme_npoints
1140 #normals = np.einsum('ij,jik->ik', hull.A, x)
1141 np.einsum('ij,jik->ik', hull.A, x, out=r) # normals
1142 result *= r
1091 1143
1092 1144 # n_values x nsimplex x scheme_points # n_values x nsimplex x scheme_points
1093 return cdfs / 2 / r / np.pi * normals
1145 return result
1094 1146
1095 1147
1148 def _chi_outside_callback(hull, x):
1149 # x.shape == (simplex_dim + 1, nsimplex, scheme_npoints) #(3, 26, 56)
1150
1151 nvar = hull.sample.nvar
1152
1153 #č přijde velké, obrovské numpy pole
1154 #č a já bych chtěl, aby se tu zbytečně nealokovalo
1155
1156 result = np.einsum('ijk,ijk->jk', x, x) # r**2
1157 r_nvar = np.power(result, nvar/2) # r**(nvar-1) (chi) * r (to normalize)
1158
1159 np.divide(result, 2, out=result) # r**2 / 2
1160 special.gammaincc(nvar / 2, result, out=result) # F
1161
1162 #č normalizujeme normály. Jenže normály samotné dodáme později
1163 result /= r_nvar
1164
1165 # A = nsimplex x nvar
1166 # x = nvar x nsimplex x scheme_npoints
1167 # normals = nsimplex x scheme_npoints
1168 #normals = np.einsum('ij,jik->ik', hull.A, x)
1169 np.einsum('ij,jik->ik', hull.A, x, out=r_nvar) # normals
1170 result *= r_nvar
1171
1172 # n_values x nsimplex x scheme_points
1173 return result
1096 1174
1097 1175 def get_convex_hull_estimation(hull): def get_convex_hull_estimation(hull):
1098 1176
Hints:
Before first commit, do not forget to setup your git environment:
git config --global user.name "your_name_here"
git config --global user.email "your@email_here"

Clone this repository using HTTP(S):
git clone https://rocketgit.com/user/iam-git/WellMet

Clone this repository using ssh (do not forget to upload a key first):
git clone ssh://rocketgit@ssh.rocketgit.com/user/iam-git/WellMet

Clone this repository using git:
git clone git://git.rocketgit.com/user/iam-git/WellMet

You are allowed to anonymously push to this repository.
This means that your pushed commits will automatically be transformed into a merge request:
... clone the repository ...
... make some changes and some commits ...
git push origin main