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)
sball: new gamma-based solution. Keep the old version as well b6ac51b6b08ca2e7d5c398fa56c5a9a0832ab301 I am 2021-02-24 10:55:11
mart: add convergence plot 85a9ef727d15fc6567b7b0a594ed52cef9fe5680 I am 2021-02-22 23:32:34
mart: add basic plot and scatter functions 20e04c27db7ee3ac3a6418300eb4fc1ef532a1c7 I am 2021-02-22 14:45:31
mart: add convex hull related functions 5284e4c1c747fe25d220e7e3a6a9803d8bc6d4a1 I am 2021-02-21 20:57:52
axes3d_: draw only bottom pane 5a24020e052a7adbc17485b4c9c74d7eea801765 Alex 2021-02-12 20:50:09
mart3d: tri surface fix 2f1889070cea6c97937e00d76ce029744c9f8d07 Alex 2021-02-12 10:01:59
mart3d: add plot_wireframe 8b21a8ca659e4f2146459ea5f4cc006535f51308 I am 2021-02-12 04:45:56
mart and mart3d: add new functions fdbfd002214bf7ab7be8df8e90d6b96b62c82d13 I am 2021-02-10 17:18:02
qt_plot: update slider on redraw, add simplex estimation with no graphics 7b392cd9ee8487a620466fe757a17c7fcfff2770 I am 2021-02-09 23:00:32
mart: WIP 81d970b1e515032f0414d1cc4c6d7473f187d162 I am 2021-02-09 19:38:36
mart3d: add scatter_points() 26644c64e3b996b68e1ee133e4531380d4e00090 I am 2021-02-08 18:03:51
mart3d: new matplotlib-support module 8b6b13937a0694998e0b4de7b8a96d89abecc64a I am 2021-02-08 13:47:56
gmodels: increased no of discretization points in boundary for Pareto function 7f12597e2d15b2564d612a014fd75492663043dd Miroslav Vořechovský 2021-02-05 15:28:32
gmodels.py added boundary for Pareto function e7bc61e9c5854bc4cff0ee05b651b4b7ede4e0a1 Miroslav Vořechovský 2021-02-05 15:14:52
improved estimation of pf for Pareto Example 427602319f6a81e896b7aa50fd8947c5a1a91873 Miroslav Vořechovský 2021-02-05 14:15:41
dicebox.Razitko: little fixes 0e3f9147a40b7008ddec0893c3965136bac7fedf I am 2021-02-05 12:32:27
f_models.UnCorD.sample_pdf: Rn pdf fix 66cf81e4fc0bb2c6dab506cef5309c5da65a7f2c I am 2021-02-05 12:31:22
dicebox.Razitko: TRI_overall_estimations workaround 4b8146b788d0a857330c1894b5cf1c8d332103ce I am 2021-02-05 09:13:16
dicebox.Razitko: new DiceBox d9c1c9cdec94601bce9e81861ce1d4dfb5393481 I am 2021-02-05 08:04:28
simplex.Triangulation: little changes dc1efc3a081073d57be0f7a06330eb97f9634447 I am 2021-02-05 08:02:57
Commit b6ac51b6b08ca2e7d5c398fa56c5a9a0832ab301 - sball: new gamma-based solution. Keep the old version as well
Author: I am
Author date (UTC): 2021-02-24 10:55
Committer name: I am
Committer date (UTC): 2021-02-24 10:55
Parent(s): 85a9ef727d15fc6567b7b0a594ed52cef9fe5680
Signer:
Signing key:
Signing status: N
Tree: 5487696f5be664d26185d5a88f5085baa76f2408
File Lines added Lines deleted
sball.py 75 124
sball_old.py 0 0
File sball.py changed (mode: 100644) (index bac994a..e9ebd21)
2 2 # coding: utf-8 # coding: utf-8
3 3
4 4 import numpy as np import numpy as np
5 #from scipy import stats
6 from scipy import special # for S_ball
7 from scipy import integrate # for S_ball
5 import scipy.special as sc
8 6
9 # нельзя просто так взять и написать Ньютонову методу
10 # fails on nvar = 501, fails on Sball(500).get_r(0), fails on Sball(800).get_r(0.999)
11 7
12 class Sball:
8 #######################################################
9 # s-balls  -- tools to calc probabilities and radii ###
10 #######################################################
11
12
13 def get_ps_ball(d,R):
14 "returns probability of falling inside d-ball with radius R"
15 #return np.sum(np.exp(-(rho**2)/2)*rho**(d-1) )* R/n
16 return sc.gammainc(d/2, R**2/2)
17
18 def get_pf_ball(d,R):
19 "returns probability of falling outside d-ball with radius R"
20 return sc.gammaincc(d/2, R**2/2)
21
22 def get_Radius_pf(d,pf_ball):
23 "returns radius of an d-ball with probability pf outside"
24 rsqdiv2 = sc.gammainccinv(d/2, pf_ball)
25 return np.sqrt(2*rsqdiv2) #radius
26
27 def get_Radius_ps(d,ps_ball):
28 "returns radius of an d-ball with probability ps inside"
29 rsqdiv2 = sc.gammaincinv(d/2, ps_ball)
30 return np.sqrt(2*rsqdiv2) #radius
31
32
33
34
35
36 # implement class compatible to the old ones
37
38 # dispatcher
39 def Sball(nvar):
40 if nvar == 2:
41 return Sball_2D(nvar)
42 else:
43 return Sball_nD(nvar)
44
45 class Sball_nD:
13 46 def __init__(self, nvar): def __init__(self, nvar):
14 47 self.nvar = nvar self.nvar = nvar
15 if nvar != 2:
16 self.C = 2**(1-nvar/2) / special.gamma(nvar/2)
17 self.logC = (1-nvar/2)*np.log(2) - special.gammaln(nvar/2)
18 self.flex = self.current_r = np.sqrt(self.nvar-1)
19 self.flex_pf = self.current_pf = self.get_pf(self.flex)
20
48 self.a = nvar/2
49
21 50 def get_pf(self, r): def get_pf(self, r):
22 """
23 returns pf, i.e. complementary part of multidimensional Gaussian distribution
24 """
25 if self.nvar == 1:
26 #return 1 - 2**(1-nvar/2) / special.gamma(nvar/2) * (np.sqrt(np.pi)*special.erf(r/np.sqrt(2)))/np.sqrt(2)
27 return 1 - special.erf(r/1.4142135623730951)
28 elif self.nvar == 2:
29 return np.exp(-r**2/2)
30 elif self.nvar == 3:
31 #return 1 - 2**(1-nvar/2) / special.gamma(nvar/2) * (np.exp(-r**2/2)*(np.sqrt(np.pi)*np.exp(r**2/2)*special.erf(r/np.sqrt(2))-np.sqrt(2)*r))/np.sqrt(2)
32 return 1 - 0.5641895835477564 * (np.exp(-r**2/2)*(np.sqrt(np.pi)*np.exp(r**2/2)*special.erf(r/np.sqrt(2))-np.sqrt(2)*r))
33 elif self.nvar == 4:
34 return (r**2/2+1)*np.exp(-r**2/2)
35 elif self.nvar == 6:
36 return (r**4+4*r**2+8)*np.exp(-r**2/2)/8
37
38 # nvar=8: (48-(r^6+6*r^4+24*r^2+48)*e^(-r^2/2) / 2**(nvar/2))/48
39
40 # hračička ve hračce
41 # nemám žádnou jistotu, že tohle počítá přesněji
42 # ale ve výsokých dimenzích aspoň počítá
43 elif self.nvar % 2 == 0: # sudé
44 poly = [1]
45 for i in range(self.nvar-2, 0, -2):
46 poly.append(0)
47 poly.append(i*poly[-2])
48 return np.polyval(np.array(poly) / poly[-1], r) * np.exp(-r**2/2)
49
50 else:
51 try:
52 pf = self.C * integrate.quad(lambda x: np.exp(-(x**2)/2)*x**(self.nvar-1), r, np.inf)[0]
53 except OverflowError:
54 pf = 1 - self.C * integrate.quad(lambda x: np.exp(-(x**2)/2)*x**(self.nvar-1), 0, r)[0]
55
56 return pf
57
51 "returns pf, i.e. complementary part of multidimensional Gaussian distribution"
52 return sc.gammaincc(self.a, r**2/2)
53
54 def get_ps(self, r):
55 "returns probability of falling inside d-ball with radius R"
56 return sc.gammainc(self.a, r**2/2)
57
58 58 def get_r(self, desired_pf): def get_r(self, desired_pf):
59 """
60 sball_inversion
61 returns r
62 """
63 if self.nvar == 2:
64 return np.sqrt(-2*np.log(desired_pf))
65 elif self.flex_pf == desired_pf:
66 return self.flex
67 else:
68 # je to jistější
69 self.current_r = self.flex
70 self.current_pf = previous_pf = self.flex_pf
71
72 self.__do_iter(desired_pf)
73 self.current_pf = self.get_pf(self.current_r)
74 # hrůza
75 # pokračujeme, dokud to nezkonverguje, přenejmenším pokud to konvergue a neosciluje.
76 while self.current_pf != previous_pf and self.current_pf != desired_pf\
77 and (self.current_pf > desired_pf or previous_pf < desired_pf):
78
79 previous_pf = self.current_pf
80 self.__do_iter(desired_pf)
81 self.current_pf = self.get_pf(self.current_r)
82
83 return self.current_r
84
85 def __do_iter(self, desired_pf):
86 r = self.current_r
87 denominator = (self.C * np.exp(-(r**2)/2)*r**(self.nvar-1))
88 if denominator != 0 and not np.isnan(denominator):
89 self.current_r += (self.current_pf - desired_pf) / denominator
90 else:
91 # zkombinujeme C a r^nvar, ale stejně nikoho to nezahraní
92 log_delta = np.log(abs(self.current_pf - desired_pf)) + (r**2)/2 - (np.log(r)*(self.nvar-1) + self.logC)
93 self.current_r += np.exp(log_delta)*np.sign(self.current_pf - desired_pf)
94
95 if self.current_r < 0:
96 self.current_r = r/2
97
98
59 "sball inversion. Returns radius of the s-ball with probability pf outside"
60 rsqdiv2 = sc.gammainccinv(self.a, desired_pf)
61 return np.sqrt(2*rsqdiv2) #radius
62
63
99 64 def get_r_iteration(self, desired_pf): def get_r_iteration(self, desired_pf):
100 """
101 Same as get_r, but do just one iteration
102 """
103
104
105 if self.nvar == 2:
106 return np.sqrt(-2*np.log(desired_pf)), desired_pf
107
108 # logaritmus je na nulu citelný
109 elif self.current_pf - desired_pf != 0:
110
111 # hrůza, nečitelný
112 # pokud je současné r-ko v jiné straně od chtěného r-ka, tak se vrátíme do inflexního bodu
113 if (self.flex_pf > self.current_pf) is (self.flex_pf < desired_pf):
114 # vstupní kontrola
115 self.current_r = self.flex
116 self.current_pf = self.flex_pf
117
118
119 r = self.current_r # pro výstupní kontrolu
120
121 self.__do_iter(desired_pf)
122
123
124 # vystupní kontrola
125 if (self.flex > self.current_r) is (self.flex < r):
126 # preskočili jsme inflexní bod
127 self.current_r = self.flex
128 self.current_pf = self.flex_pf
129 # ještě jednou
130 self.__do_iter(desired_pf)
131
132 self.current_pf = self.get_pf(self.current_r) # ne že bychom pf potrebovali v tomto kroce, ale...
133 return self.current_r, self.current_pf
134
135
65 "Same as .get_r(), just keeps compatibility with previous versions"
66 return self.get_r(desired_pf), desired_pf
67
68
69
70 class Sball_2D(Sball_nD):
71 def get_pf(self, r):
72 "returns pf, i.e. complementary part of multidimensional Gaussian distribution"
73 return np.exp(-r**2/2)
74
75 def get_r(self, desired_pf):
76 "sball inversion. Returns radius of the s-ball with probability pf outside"
77 return np.sqrt(-2*np.log(desired_pf))
78
79
80 # calculation is as fast as Sball_nD
81 # but I'm not sure about precision
82 class Sball_4D(Sball_nD):
83 def get_pf(self, r):
84 "returns pf, i.e. complementary part of multidimensional Gaussian distribution"
85 return (r**2/2+1)*np.exp(-r**2/2)
86
136 87
137 88
File sball_old.py copied from file sball.py (similarity 100%)
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