From 8b4c0c6dd4d8612ef2614f73be1e55d82a19227b Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Tue, 6 Jun 2017 05:01:40 +0200 Subject: [PATCH] Added markov.py --- wafo/markov.py | 1431 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1431 insertions(+) create mode 100644 wafo/markov.py diff --git a/wafo/markov.py b/wafo/markov.py new file mode 100644 index 0000000..ced5cc2 --- /dev/null +++ b/wafo/markov.py @@ -0,0 +1,1431 @@ +from __future__ import division +import numpy as np + + +def nt2fr(nt, kind=0): + """ + NT2FR Calculates the frequency matrix given the counting distribution matrix. + + Parameters + ---------- + nt = the counting distribution matrix, + kind = 0,1 + + Returns + ------- + fr = the frequency matrix, + + If kind=0 function computes the inverse to + + N_T(u,v) = #{ (M_i,m_i); M_i>u, m_i=u, m_i= epsilon + if converged: + break + + F = np.fliplr(fmM) + frfc = np.fliplr(fun2rfc(fmM)) + + return F, frfc + + +def iter_mc(frfc, fmM_0=None, k=1, epsilon=1e-5): + """ + ITER_MC Calculates a kernel of a MC given a rainflow matrix + + Solves f_rfc = f_xy + F_mc(f_xy) for f_xy. + + Call: [fmM_k frfc_k]=iter_mc(frfc,fmM_0,k,eps) + + fmM_k = the solution to the equation frfc = fmM + F(fmM), + frfc_k = the rainflow matrix; frfc_k = fmM_k + F(fmM_k). + + + frfc = the rainflow matrix to be inverted, + fmM_0 = the first approximation to the Markov matrix, if not + specified fmM_0=frfc, + k = number of iterations, if not specified, k=1. + eps = a convergence treshold, default value; eps=0.00001 + + See also + -------- + iter_, spec2cmat, mctp2rfm, mc2rfm + + References + ---------- + Rychlik, I. (1996) + 'Simulation of load sequences from Rainflow matrices: Markov method' + Int. J. Fatigue, Vol 18, pp 429-438 + """ + return _raw_iter(mc2rfc, frfc, fmM_0, k, epsilon) + + +def _raise_kind_error(kind): + if kind in (-1, 0): + raise NotImplementedError('kind = {} not yet implemented'.format(kind)) + else: + raise ValueError('kind = {}: not a valid value of kind'.format(kind)) + + +def nt2cmat(nt, kind=1): + """ + Return cycle matrix from a counting distribution. + + Parameters + ---------- + NT: 2D array + Counting distribution. [nxn] + kind = 1: causes peaks to be projected upwards and troughs + downwards to the closest discrete level (default). + = 0: causes peaks and troughs to be projected to + the closest discrete level. + = -1: causes peaks to be projected downwards and the + troughs upwards to the closest discrete level. + + Returns + ------- + cmat = Cycle matrix. [nxn] + + Example + -------- + >>> import numpy as np + >>> cmat0 = np.round(np.triu(np.random.rand(4, 4), 1)*10) + >>> cmat0 = np.array([[ 0., 5., 6., 9.], + ... [ 0., 0., 1., 7.], + ... [ 0., 0., 0., 4.], + ... [ 0., 0., 0., 0.]]) + + >>> nt = cmat2nt(cmat0) + >>> np.allclose(nt, + ... [[ 0., 0., 0., 0.], + ... [ 20., 15., 9., 0.], + ... [ 28., 23., 16., 0.], + ... [ 32., 27., 20., 0.]]) + True + >>> cmat = nt2cmat(nt) + >>> np.allclose(cmat, [[ 0., 5., 6., 9.], + ... [ 0., 0., 1., 7.], + ... [ 0., 0., 0., 4.], + ... [ 0., 0., 0., 0.]]) + True + + See also + -------- + cmat2nt + """ + n = len(nt) # Number of discrete levels + if kind == 1: + I = np.r_[0:n - 1] + J = np.r_[1:n] + c = nt[I+1][:, J-1] - nt[I][:, J-1] - nt[I+1][:, J] + nt[I][:, J] + c2 = np.vstack((c, np.zeros((n - 1)))) + cmat = np.hstack((np.zeros((n, 1)), c2)) + elif kind == 11: # same as def=1 but using for-loop + cmat = np.zeros((n, n)) + j = np.r_[1:n] + for i in range(n - 1): + cmat[i, j] = nt[i+1, j-1] - nt[i, j-1] - nt[i+1, j] + nt[i, j] + else: + _raise_kind_error(kind) + return cmat + + +def cmat2nt(cmat, kind=1): + """ + CMAT2NT Calculates a counting distribution from a cycle matrix. + + Parameters + ---------- + cmat = Cycle matrix. [nxn] + kind = 1: causes peaks to be projected upwards and troughs + downwards to the closest discrete level (default). + = 0: causes peaks and troughs to be projected to + the closest discrete level. + = -1: causes peaks to be projected downwards and the + troughs upwards to the closest discrete level. + Returns + ------- + NT: n x n array + Counting distribution. + + Example + ------- + >>> import numpy as np + >>> cmat0 = np.round(np.triu(np.random.rand(4, 4), 1)*10) + >>> cmat0 = np.array([[ 0., 5., 6., 9.], + ... [ 0., 0., 1., 7.], + ... [ 0., 0., 0., 4.], + ... [ 0., 0., 0., 0.]]) + + >>> nt = cmat2nt(cmat0, kind=11) + >>> np.allclose(nt, + ... [[ 0., 0., 0., 0.], + ... [ 20., 15., 9., 0.], + ... [ 28., 23., 16., 0.], + ... [ 32., 27., 20., 0.]]) + True + + >>> cmat = nt2cmat(nt, kind=11) + >>> np.allclose(cmat, [[ 0., 5., 6., 9.], + ... [ 0., 0., 1., 7.], + ... [ 0., 0., 0., 4.], + ... [ 0., 0., 0., 0.]]) + True + + See also + -------- + nt2cmat + """ + n = len(cmat) # Number of discrete levels + nt = np.zeros((n, n)) + + if kind == 1: + csum = np.cumsum + flip = np.fliplr + nt[1:n, :n - 1] = flip(csum(flip(csum(cmat[:-1, 1:], axis=0)), axis=1)) + elif kind == 11: # same as def=1 but using for-loop + # j = np.r_[1:n] + for i in range(1, n): + for j in range(n - 1): + nt[i, j] = np.sum(cmat[:i, j + 1:n]) + else: + _raise_kind_error(kind) + return nt + + +def mctp2tc(f_Mm, utc, param, f_mM=None): + """ + MCTP2TC Calculates frequencies for the upcrossing troughs and crests + using Markov chain of turning points. + + Parameters + ---------- + f_Mm = the frequency matrix for the Max2min cycles, + utc = the reference level, + param = a vector defining the discretization used to compute f_Mm, + note that f_mM has to be computed on the same grid as f_mM. + f_mM = the frequency matrix for the min2Max cycles. + + Returns + ------- + f_tc = the matrix with frequences of upcrossing troughs and crests, + + Example + ------- + >>> fmM = np.array([[ 0.0183, 0.0160, 0.0002, 0.0000, 0], + ... [0.0178, 0.5405, 0.0952, 0, 0], + ... [0.0002, 0.0813, 0, 0, 0], + ... [0.0000, 0, 0, 0, 0], + ... [ 0, 0, 0, 0, 0]]) + >>> param = (-1, 1, len(fmM)) + >>> utc = 0 + >>> f_tc = mctp2tc(fmM, utc, param) + >>> np.allclose(f_tc, + ... [[ 0.0, 1.59878359e-02, -1.87345256e-04, 0.0, 0.0], + ... [ 0.0, 5.40312726e-01, 3.86782958e-04, 0.0, 0.0], + ... [ 0.0, 0.0, 0.0, 0.0, 0.0], + ... [ 0.0, 0.0, 0.0, 0.0, 0.0], + ... [ 0.0, 0.0, 0.0, 0.0, 0.0]]) + True + """ + def _check_ntc(ntc, n): + if ntc > n - 1: + raise IndexError('index for mean-level out of range, stop') + + def _check_discretization(param, ntc): + if not (1 < ntc < param[2]): + raise ValueError('the reference level out of range, stop') + + def _normalize_rows(arr): + n = len(arr) + for i in range(n): + rowsum = np.sum(arr[i]) + if rowsum != 0: + arr[i] = arr[i] / rowsum + return arr + + def _make_tempp(P, Ph, i, ntc): + Ap = P[i:ntc - 1, i + 1:ntc] + Bp = Ph[i + 1:ntc, i:ntc - 1] + dim_p = ntc - 1 - i + tempp = np.zeros((dim_p, 1)) + I = np.eye(len(Ap)) + if i == 1: + e = Ph[i + 1:ntc, 0] + else: + e = np.sum(Ph[i + 1:ntc, :i - 1], axis=1) + + if max(abs(e)) > 1e-10: + if dim_p == 1: + tempp[0] = (Ap / (1 - Bp * Ap) * e) + else: + rh = I - np.dot(Bp, Ap) + tempp = np.dot(Ap, np.linalg.solve(rh, e)) + + # end + # end + return tempp + + def _make_tempm(P, Ph, j, ntc, n): + Am = P[ntc-1:j, ntc:j+1] + Bm = Ph[ntc:j+1, ntc-1:j] + dim_m = j - ntc + 1 + tempm = np.zeros((dim_m, 1)) + Im = np.eye(len(Am)) + if j == n - 1: + em = P[ntc-1:j, n] + else: + em = np.sum(P[ntc-1:j, j + 1:n], axis=1) + # end + if max(abs(em)) > 1e-10: + if dim_m == 1: + tempm[0] = (Bm / (1 - Am * Bm) * em) + else: + rh = Im - np.dot(Am, Bm) + tempm = np.dot(Bm, np.linalg.lstsq(rh, em)[0]) + # end + # end + return tempm + + if f_mM is None: + f_mM = np.copy(f_Mm) * 1.0 + + u = np.linspace(*param) + udisc = u[::-1] # np.fliplr(u) + ntc = np.sum(udisc >= utc) + n = len(f_Mm) + _check_ntc(ntc, n) + _check_discretization(param, ntc) + + # normalization of frequency matrices + f_Mm = _normalize_rows(f_Mm) + P = np.fliplr(f_Mm) + Ph = np.rot90(np.fliplr(f_mM*1.0), -1) + Ph = _normalize_rows(Ph) + Ph = np.fliplr(Ph) + + F = np.zeros((n, n)) + F[:ntc - 1, :(n - ntc)] = f_mM[:ntc - 1, :(n - ntc)] + F = cmat2nt(F) + + for i in range(1, ntc): + for j in range(ntc-1, n - 1): + if i < ntc-1: + tempp = _make_tempp(P, Ph, i, ntc) + b = np.dot(np.dot(tempp.T, f_mM[i:ntc - 1, n - j - 2::-1]), + np.ones((n - j - 1, 1))) + # end + if j > ntc-1: + tempm = _make_tempm(P, Ph, j, ntc, n) + c = np.dot(np.dot(np.ones((1, i)), + f_mM[:i, n - ntc-1:n - j - 2:-1]), + tempm) + # end + if (j > ntc-1) and (i < ntc-1): + a = np.dot(np.dot(tempp.T, + f_mM[i:ntc - 1, n - ntc - 1:-1:n - j + 1]), + tempm) + F[i, n - j - 1] = F[i, n - j - 1] + a + b + c + # end + if (j == ntc-1) and (i < ntc-1): + F[i, n - ntc] = F[i, n - ntc] + b + for k in range(ntc): + F[i, n - k - 1] = F[i, n - ntc] + # end + # end + if (j > ntc-1) and (i == ntc-1): + F[i, n - j - 1] = F[i, n - j - 1] + c + for k in range(ntc-1, n): + F[k, n - j - 1] = F[ntc-1, n - j - 1] + # end + # end + # end + # end + + # fmax=max(max(F)); + # contour (u,u,flipud(F),... + # fmax*[0.005 0.01 0.02 0.05 0.1 0.2 0.4 0.6 0.8]) + # axis([param(1) param(2) param(1) param(2)]) + # title('Crest-trough density') + # ylabel('crest'), xlabel('trough') + # axis('square') + # if mlver>1, commers, end + return nt2cmat(F) + + +def mctp2rfc(fmM, fMm=None): + ''' + Return Rainflow matrix given a Markov chain of turning points + + computes f_rfc = f_mM + F_mct(f_mM). + + Parameters + ---------- + fmM = the min2max Markov matrix, + fMm = the max2min Markov matrix, + + Returns + ------- + f_rfc = the rainflow matrix, + + Example: + ------- + >>> fmM = np.array([[ 0.0183, 0.0160, 0.0002, 0.0000, 0], + ... [0.0178, 0.5405, 0.0952, 0, 0], + ... [0.0002, 0.0813, 0, 0, 0], + ... [0.0000, 0, 0, 0, 0], + ... [ 0, 0, 0, 0, 0]]) + + >>> np.allclose(mctp2rfc(fmM), + ... [[ 2.669981e-02, 7.799700e-03, 4.906077e-07, 0.0, 0.0], + ... [ 9.599629e-03, 5.485009e-01, 9.539951e-02, 0.0, 0.0], + ... [ 5.622974e-07, 8.149944e-02, 0.0, 0.0, 0.0], + ... [ 0.0, 0.0, 0.0, 0.0, 0.0], + ... [ 0.0, 0.0, 0.0, 0.0, 0.0]], 1.e-7) + True + ''' + def _get_PMm(AA1, MA, nA): + PMm = AA1.copy() + for j in range(nA): + norm = MA[j] + if norm != 0: + PMm[j, :] = PMm[j, :] / norm + PMm = np.fliplr(PMm) + return PMm + + if fMm is None: + fmM = np.atleast_1d(fmM) + fMm = fmM + else: + fmM, fMm = np.atleast_1d(fmM, fMm) + f_mM, f_Mm = fmM.copy(), fMm.copy() + N = max(f_mM.shape) + f_max = np.sum(f_mM, axis=1) + f_min = np.sum(f_mM, axis=0) + f_rfc = np.zeros((N, N)) + f_rfc[N - 2, 0] = f_max[N - 2] + f_rfc[0, N - 2] = f_min[N - 2] + for k in range(2, N - 1): + for i in range(1, k): + AA = f_mM[N - 1 - k:N - 1 - k + i, k - i:k] + AA1 = f_Mm[N - 1 - k:N - 1 - k + i, k - i:k] + RAA = f_rfc[N - 1 - k:N - 1 - k + i, k - i:k] + nA = max(AA.shape) + MA = f_max[N - 1 - k:N - 1 - k + i] + mA = f_min[k - i:k] + SA = AA.sum() + SRA = RAA.sum() + + DRFC = SA - SRA + NT = min(mA[0] - sum(RAA[:, 0]), MA[0] - sum(RAA[0, :])) # check! + NT = max(NT, 0) # ??check + + if NT > 1e-6 * max(MA[0], mA[0]): + NN = MA - np.sum(AA, axis=1) # T + e = (mA - np.sum(AA, axis=0)) # T + e = np.flipud(e) + PmM = np.rot90(AA.copy()) + for j in range(nA): + norm = mA[nA - 1 - j] + if norm != 0: + PmM[j, :] = PmM[j, :] / norm + e[j] = e[j] / norm + # end + # end + fx = 0.0 + if (max(np.abs(e)) > 1e-6 and + max(np.abs(NN)) > 1e-6 * max(MA[0], mA[0])): + PMm = _get_PMm(AA1, MA, nA) + + A = PMm + B = PmM + + if nA == 1: + fx = NN * (A / (1 - B * A) * e) + else: + rh = np.eye(A.shape[0]) - np.dot(B, A) + # least squares + fx = np.dot(NN, np.dot(A, np.linalg.solve(rh, e))) + # end + # end + f_rfc[N - 1 - k, k - i] = fx + DRFC + + # check2=[ DRFC fx] + # pause + else: + f_rfc[N - 1 - k, k - i] = 0.0 + # end + # end + m0 = max(0, f_min[0] - np.sum(f_rfc[N - k + 1:N, 0])) + M0 = max(0, f_max[N - 1 - k] - np.sum(f_rfc[N - 1 - k, 1:k])) + f_rfc[N - 1 - k, 0] = min(m0, M0) + # n_loops_left=N-k+1 + # end + + for k in range(1, N): + M0 = max(0, f_max[0] - np.sum(f_rfc[0, N - k:N])) + m0 = max(0, f_min[N - 1 - k] - np.sum(f_rfc[1:k + 1, N - 1 - k])) + f_rfc[0, N - 1 - k] = min(m0, M0) + # end + +# clf +# subplot(1,2,2) +# pcolor(levels(paramm),levels(paramM),flipud(f_mM)) +# title('Markov matrix') +# ylabel('max'), xlabel('min') +# axis([paramm(1) paramm(2) paramM(1) paramM(2)]) +# axis('square') +# +# subplot(1,2,1) +# pcolor(levels(paramm),levels(paramM),flipud(f_rfc)) +# title('Rainflow matrix') +# ylabel('max'), xlabel('rfc-min') +# axis([paramm(1) paramm(2) paramM(1) paramM(2)]) +# axis('square') + + return f_rfc + + +def mc2rfc(f_xy, paramv=None, paramu=None): + """ + MC2RFC Calculates a rainflow matrix given a Markov chain with kernel f_xy; + f_rfc = f_xy + F_mc(f_xy). + + CALL: f_rfc = mc2rfc(f_xy); + + where + + f_rfc = the rainflow matrix, + f_xy = the frequency matrix of Markov chain (X0,X1) + but only the triangular part for X1>X0. + + Further optional input arguments; + + CALL: f_rfc = mc2rfc(f_xy,paramx,paramy); + + paramx = the parameter matrix defining discretization of x-values, + paramy = the parameter matrix defining discretization of y-values, + """ + N = len(f_xy) + if paramv is None: + paramv = (-1, 1, N) + + + if paramu is None: + paramu = paramv + + dd = np.diag(np.rot90(f_xy)) + Splus = np.sum(f_xy, axis=1).T + Sminus = np.fliplr(sum(f_xy)) + Max_rfc = np.zeros((N, 1)) + Min_rfc = np.zeros((N, 1)) + norm = np.zeros((N, 1)) + for i in range(N): + Spm=Sminus[i] + Splus[i]-dd[i] + if Spm > 0: + Max_rfc[i]=(Splus[i]-dd[i])*(Splus[i]-dd[i])/(1-dd[i]/Spm)/Spm + Min_rfc[i]=(Sminus[i]-dd[i])*(Sminus[i]-dd[i])/(1-dd[i]/Spm)/Spm + norm[i]=Spm + # end if + # end for + + #cross=zeros(N,1) + #for i=2:N + # cross(N-i+1)=cross(N-i+2)+Sminus(N-i+2)-Splus(N-i+2) + #end + + f_rfc= np.zeros((N,N)) + f_rfc[N-1, 1] = Max_rfc[N-1] + f_rfc[1, N-1] = Min_rfc[2] + + for k in range(2, N-1): + for i in range(1, k-1): + + # AAe= f_xy(1:N-k,1:k-i) + # SAAe=sum(sum(AAe)) + AA = f_xy[N-k:N-k+i, :][:, k-i:k] + #RAA = f_rfc(N-k+1:N-k+i,k-i+1:k) + RAA = f_rfc[N-k:N-k+i, :][:, k-i:k] + nA = len(AA) + # MA = Splus(N-k+1:N-k+i) + MA = Splus[N-k:N-k+i] + mA = Sminus[N-k:N-k+i] + normA = norm[N-k:N-k+i] + MA_rfc = Max_rfc[N-k:N-k+i] + # mA_rfc=Min_rfc(k-i+1:k) + SA = np.sum(AA) + SRA = np.sum(RAA) + SMA_rfc = sum(MA_rfc) + SMA = np.sum(MA) + DRFC = SA-SMA-SRA+SMA_rfc + + NT = MA_rfc[0] - np.sum(RAA[0, :]) + + # if k==35 + # check=[MA_rfc(1) sum(RAA(1,:))] + # pause + # end + + NT = np.maximum(NT, 0) + + if NT > 1e-6*MA_rfc[0]: + + NN = MA - np.sum(AA,axis=1).T + e = (np.fliplr(mA)-np.sum(AA)).T + e = np.flipud(e) + AA = AA + np.flipud(np.rot90(AA,-1)) + AA = np.rot90(AA) + AA = AA - 0.5*np.diag(np.diag(AA)) + + + for j in range(nA): + if normA[j] !=0: + AA[j,:] = AA[j, :] / normA[j] + e[j] = e[j] / normA[j] + # end if + # end for + fx = 0. + + if np.max(np.abs(e)) > 1e-7 and np.max(np.abs(NN)) > 1e-7*MA_rfc[0]: + I = np.eye(np.shape(AA)) + + if nA==1: + fx = NN/(1-AA)*e + else: + # TODO CHECK this + fx = NN*np.linalg.solve((I-AA), e)[0] # (I-AA)\e + # end + # end + + f_rfc[N-k, k-i] = DRFC + fx + #end + # end + m0 = np.maximum(0, Min_rfc[N]-sum(f_rfc[N-k+1:N, 0])) + M0 = np.maximum(0, Max_rfc[N-k]-sum(f_rfc[N-k, 1:k])); + f_rfc[N-k, 0] = min(m0, M0) + # n_loops_left=N-k+1 + # end for + + for k in range(1, N): + M0 = max(0, Max_rfc[0] -sum(f_rfc[0, N-k+1:N])) + m0 = max(0, Min_rfc[k] - sum(f_rfc[1:k, N-k])) + f_rfc[0, N-k] = min(m0, M0) + # end for + f_rfc = f_rfc + np.rot90(np.diag(dd),-1) + +# clf +# subplot(1,2,2) +# pcolor(levels(paramv),levels(paramu),flipud(f_xy+flipud(rot90(f_xy,-1)))) +# axis([paramv(1), paramv(2), paramu(1), paramu(2)]) +# title('MC-kernel f(x,y)') +# ylabel('y'), xlabel('x') +# axis('square') +# +# subplot(1,2,1) +# pcolor(levels(paramv),levels(paramu),flipud(f_rfc)) +# axis([paramv(1), paramv(2), paramu(1), paramu(2)]) +# title('Rainflow matrix') +# ylabel('max'), xlabel('rfc-min') +# axis('square') + + return f_rfc + + +def mktestmat(param=(-1, 1, 32), x0=None, s=None, lam=1, numsubzero=0): + """ + + MKTESTMAT Makes test matrices for min-max (and max-min) matrices. + + Parameters + ---------- + param = Parameter vector, [a b n], defines discretization. + x0 = Center of ellipse. [min Max] [1x2] + s = Standard deviation. (0>> import numpy as np + >>> param = [-1, 1, 32] + >>> F, Fh = mktestmat(param, x0=[-0.2, 0.2], s=0.25, lam=0.5) + >>> np.allclose(F[0,:4], + ... [ 8.45677332e-29, 5.94193456e-27, 3.53462989e-25, 1.78013497e-23]) + True + + >>> u = np.linspace(*param) + + cmatplot(u,u,F,3) + axis('square') + >>> F, Fh = mktestmat(param, x0=[-0.2, 0.2], s=0.25, lam=0.5, numsubzero=-np.inf) + >>> np.allclose(F[0,:4], + ... [ 8.45677332e-29, 5.94193456e-27, 3.53462989e-25, 1.78013497e-23]) + True + + cmatplot(u,u,F,3); axis('square'); + + close all; + """ +# History: +# Revised by PJ 23-Nov-1999 +# updated for WAFO +# Created by PJ (Paer Johannesson) 1997 +# Copyright (c) 1997 by Paer Johannesson +# Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997 + + if x0 is None: + x0 = np.ones((2,)) * (param[1]+param[0])/2 + if s is None: + s = (param[1]-param[0])/4 + + if np.isinf(numsubzero): + numsubzero = -(param[2]+1) + + u = np.linspace(*param) + n = param[2] + + # F - min-Max matrix + + F = np.zeros((n,n)) + S = 1/2*s**2*np.array([[lam**2+1, lam**2-1], + [lam**2-1, lam**2+1]]) + Sm1 = np.linalg.pinv(S) + for i in range(min(n-1-numsubzero, n)): + for j in range(max(i + numsubzero, 0), n): + dx = np.array([u[i], u[j]]) - x0 + F[i, j] = np.exp(-1/2*np.dot(np.dot(dx, Sm1), dx)) + + Fh = F.T # Time-reversible Max-min matrix + return F, Fh + + +def cmatplot(cmat, ux=None, uy=None, method=1, clevels=None): + """ + CMATPLOT Plots a cycle matrix, e.g. a rainflow matrix. + + CALL: cmatplot(F) + cmatplot(F,method) + cmatplot(ux,uy,F) + cmatplot(ux,uy,F,method) + + F = Cycle matrix (e.g. rainflow matrix) [nxm] + method = 1: mesh-plot (default) + 2: surf-plot + 3: pcolor-plot [axis('square')] + 4: contour-plot [axis('square')] + 5: TechMath-plot [axis('square')] + 11: From-To, mesh-plot + 12: From-To, surf-plot + 13: From-To, pcolor-plot [axis('square')] + 14: From-To, contour-plot [axis('square')] + 15: From-To, TechMath-plot [axis('square')] + ux = x-axis (default: 1:m) + uy = y-axis (default: 1:n) + + Examples: + param = [-1 1 64]; u=levels(param); + F = mktestmat(param,[-0.2 0.2],0.25,1/2); + cmatplot(F,method=1); + cmatplot(u,u,F,method=2); colorbar; + cmatplot(u,u,F,method=3); colorbar; + cmatplot(u,u,F,method=4); + + close all; + + See also + -------- + cocc, plotcc + """ + + F = cmat + shape = np.shape(F) + if ux is None: + ux = np.arange(shape[1]); # Antalet kolumner + + if uy is None: + uy = np.arange(shape[0]) # Antalet rader + + if clevels is None: + Fmax = np.max(F) + if method in [5, 15]: + clevels = Fmax * np.r_[0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0] + else: # 4, 14 + clevels = Fmax*np.r_[0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8] + + # Make sure ux and uy are row vectors + ux = ux.ravel(); + uy = uy.ravel(); + + n = len(F) + + from matplotlib import pyplot as plt + if method == 1: # mesh + F = np.flipud(F.T); # Vrid cykelmatrisen for att plotta rett + plt.mesh(ux, np.fliplr(uy), F) + plt.xlabel('min') + plt.ylabel('Max') + # view(-37.5-90,30); + #v = axis; + #plt.axis([min(ux) max(ux) min(uy) max(uy) v[5:6]]); + elif method == 2: # surf + F = np.flipud(F.T) # Vrid cykelmatrisen for att plotta rett + plt.surf(ux, np.fliplr(uy),F) + plt.xlabel('min') + plt.ylabel('Max') + # view(-37.5-90,30); + #v = axis; + #plt.axis([min(ux) max(ux) min(uy) max(uy) v(5:6)]); +# elseif method == 3 # pcolor +# F = flipud(F'); +# F1 = [F zeros(length(uy),1); zeros(1,length(ux)+1)]; +# F2 = F1; F2(F2==0)=NaN; +# F1 = F2; +# dx=ux(2)-ux(1); dy=uy(2)-uy(1); +# ux1 = [ux ux(length(ux))+dx] - dx/2; +# uy1 = [uy uy(length(uy))+dy] - dy/2; +# pcolor(ux1,fliplr(uy1),F1); +# xlabel('min'); +# ylabel('Max'); +# v = axis; axis([min(ux1) max(ux1) min(uy1) max(uy1)]); +# axis('square') +# elseif method == 4 # contour +# F = flipud(F'); +# if isempty(clevels) +# Fmax=max(max(F)); +# clevels=Fmax*[0.005 0.01 0.02 0.05 0.1 0.2 0.4 0.6 0.8]; +# end +# contour(ux,fliplr(uy),F,clevels); +# xlabel('min'); +# ylabel('Max'); +# v = axis; axis([min(ux) max(ux) min(uy) max(uy)]); +# axis('square'); +# +# # Cstr=num2str(clevels(1),4); +# # for i=2:length(clevels) +# # Cstr=[Cstr ',' num2str(clevels(i),4)]; +# # end +# # title(['ISO-lines: ' Cstr]) +# +# if 1==2 +# clevels=sort(clevels); +# n_clevels=length(clevels); +# if n_clevels>12 +# disp(' Only the first 12 levels will be listed in table.') +# n_clevels=12; +# end +# +# textstart_x=0.65; +# textstart_y=0.45; +# delta_y=1/33; +# h=figtext(textstart_x,textstart_y,'Level curves at:','normalized'); +# set(h,'FontWeight','Bold') +# +# textstart_y=textstart_y-delta_y; +# +# for i=1:n_clevels +# textstart_y=textstart_y-delta_y; +# figtext(textstart_x,textstart_y,num2str(clevels(i)),'normalized') +# end +# end # 1==2 +# +# elseif method == 5 | method == 15 # TechMath-typ +# +# if isempty(clevels) +# Fmax=max(max(F)); +# clevels=Fmax*[0.001 0.005 0.01 0.05 0.1 0.5 1.0]; +# end +# v=clevels; +# # axis('ij'); +# sym = '...x+***'; +# sz = [6 20 24 8 8 8 12 16]; +# +# # plot(-1,-1,sym(1),'markersize',1),hold on +# for i = 1:length(v) +# plot(-1,-1,sym(i),'markersize',sz(i));hold on; +# end +# +# for i = 1:length(v)-1 +# Ind = (F>v(i)) & (F<=v(i+1)); +# [I,J] = find(Ind); +# # axis('ij'); +# plot(I,J,sym(i),'markersize',sz(i));hold on; +# end +# plot([1 n],[1 n],'--'); grid; +# hold off; +# +# axis([0.5 n+0.5 0.5 n+0.5]); +# +# #legendText = sprintf('#6g < f <= %6g\n',[v(1:nv-1); v(2:nv)]) +# #legendText = sprintf('<= %g\n',v(2:end)) +# +# legendText=num2str(v(1:end)'); +# +# legend(legendText,-1); +# +# title('From-To plot'); +# xlabel('To / Standing'); +# ylabel('From / Hanging'); +# +# if method == 15 +# axis('ij'); +# end +# +# elseif method == 11 # mesh +# +# mesh(ux,uy,F); +# axis('ij'); +# xlabel('To'); +# ylabel('From'); +# view(-37.5-90,30); +# v = axis; axis([min(ux) max(ux) min(uy) max(uy) v(5:6)]); +# +# elseif method == 12 # surf +# +# surf(ux,uy,F); +# axis('ij'); +# xlabel('To'); +# ylabel('From'); +# view(-37.5-90,30); +# v = axis; axis([min(ux) max(ux) min(uy) max(uy) v(5:6)]); +# +# elseif method == 13 # From-To-Matrix - pcolor +# +# F1 = [F zeros(length(uy),1); zeros(1,length(ux)+1)]; +# F2 = F1; F2(F2==0)=NaN; +# F1 = F2; +# dx=ux(2)-ux(1); dy=uy(2)-uy(1); +# ux1 = [ux ux(length(ux))+dx] - dx/2; +# uy1 = [uy uy(length(uy))+dy] - dy/2; +# axis('ij'); +# pcolor(ux1,uy1,F1) +# axis('ij'); +# xlabel('To'); +# ylabel('From'); +# v = axis; axis([min(ux1) max(ux1) min(uy1) max(uy1)]); +# axis('square') +# +# elseif method == 14 # contour +# if isempty(clevels) +# Fmax=max(max(F)); +# clevels=Fmax*[0.005 0.01 0.02 0.05 0.1 0.2 0.4 0.6 0.8]; +# end +# contour(ux,uy,F,clevels); +# axis('ij'); +# xlabel('To'); +# ylabel('From'); +# v = axis; axis([min(ux) max(ux) min(uy) max(uy)]); +# axis('square'); + +# # Cstr=num2str(clevels(1),4); +# # for i=2:length(clevels) +# # Cstr=[Cstr ',' num2str(clevels(i),4)]; +# # end +# # title(['ISO-lines: ' Cstr]) +# +# if 1==1 +# clevels=sort(clevels); +# n_clevels=length(clevels); +# if n_clevels>12 +# disp(' Only the first 12 levels will be listed in table.'); +# n_clevels=12; +# end +# +# textstart_x=0.10; +# textstart_y=0.45; +# delta_y=1/33; +# h=figtext(textstart_x,textstart_y,'Level curves at:','normalized'); +# set(h,'FontWeight','Bold') +# +# textstart_y=textstart_y-delta_y; +# +# for i=1:n_clevels +# textstart_y=textstart_y-delta_y; +# figtext(textstart_x,textstart_y,num2str(clevels(i)),'normalized') +# end +# end +# +# elseif method == 15 # TechMath-typ +# # See: 'method == 5' +# +# # if isempty(clevels) +# # Fmax=max(max(F)); +# # clevels=Fmax*[0.005 0.01 0.05 0.1 0.4 0.8]; +# # end +# # v=clevels; +# # axis('ij'); +# # sym = '...***'; +# # sz = [8 12 16 8 12 16] +# # for i = 1:length(v)-1 +# # Ind = (F>v(i)) & (F<=v(i+1)); +# # [I,J] = find(Ind); +# # axis('ij'); +# # plot(J,I,sym(i),'markersize',sz(i)),hold on +# # end +# # hold off +# # +# # axis([0.5 n+0.5 0.5 n+0.5]) +# # +# # %legendText = sprintf('%6g < f <= %6g\n',[v(1:nv-1); v(2:nv)]) +# # legendText = sprintf('<= %g\n',v(2:nv)) +# # +# # %legendText=num2str(v(2:nv)') +# # +# # legend(legendText,-1) +# +# end +# +# end % _cmatplot + +def arfm2mctp(Frfc): + """ + ARFM2MCTP Calculates the markov matrix given an asymmetric rainflow matrix. + + CALL: F = arfm2mctp(Frfc); + + F = Markov matrix (from-to-matrix) [n,n] + Frfc = Rainflow Matrix [n,n] + + Examples + -------- + param = [-1 1 32]; + u = levels(param); + F = mktestmat(param,[-0.2 0.2],0.15,2); + F = F/sum(sum(F)); + Farfc = mctp2arfm({F []}); + F1 = arfm2mctp(Farfc); + cmatplot(u,u,{F+F' F1},3); + assert(F1(20,21:25), [0.00209800691364310, 0.00266223402503216,... + 0.00300934711658560, 0.00303029619424592,... + 0.00271822008031848], 1e-10); + assert(sum(sum(abs(F1-(F+F')))), 0, 1e-10) should be zero + + close all; + + See also + -------- + rfm2mctp, mctp2arfm, smctp2arfm, cmatplot + + References + ---------- + P. Johannesson (1999): + Rainflow Analysis of Switching Markov Loads. + PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, + Lund Institute of Technology. + """ + # Tested on Matlab 5.3 + # + # History: + # Revised by PJ 09-Apr-2001 + # updated for WAFO + + # Copyright (c) 1997-1998 by Pear Johannesson + # Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998 + + # Recursive formulation a'la Igor + # + # This program used the formulation where the probabilities + # of the events are calculated using "elementary" events for + # the MCTP. + # + # Standing + # pS = Max*pS1*pS2*pS3; + # F_rfc(i,j) = pS; + # Hanging + # pH = Min*pH1*pH2*pH3; + # F_rfc(j,i) = pH; + # + # The cond. prob. pS1, pS2, pS3, pH1, pH2, pH3 are calculated using + # the elementary cond. prob. C, E, R, D, E3, Ch, Eh, Rh, Dh, E3h. + + + # T(1,:)=clock; + + N = np.sum(Frfc) + Frfc = Frfc/N; + + n = len(Frfc) # Number of levels + + # T(7,:)=clock; + # Transition matrices for MC + + Q = np.zeros((n, n)) + Qh = np.zeros((n, n)) + + # Transition matrices for time-reversed MC + + Qr = np.zeros((n, n)) + Qrh = np.zeros((n, n)) + + # Probability of minimum and of maximun + + MIN = np.sum(np.triu(Frfc).T, axis=0) + np.sum(np.tril(Frfc), axis=0) + MAX = np.sum(np.triu(Frfc), axis=0) + np.sum(np.tril(Frfc).T, axis=0) + + # Calculate rainflow matrix + + F = np.zeros((n,n)) + EYE = np.eye((n,n)) + + #fprintf(1,'Calculating row '); + for k in range(n-1): # k = subdiagonal + # fprintf(1,'-%1d',i); + + for i in range(n-k): # i = minimum + + j = i + k + 1 # maximum; + +# pS = Frfc(i,j); # Standing cycle +# pH = Frfc(j,i); # Hanging cycle +# +# Min = MIN[i] +# Max = MAX[j] +# +# # fprintf(1,'Min=%f, Max=%f\n',Min,Max); +# +# +# if j - i == 2: # Second subdiagonal +# +# # For Part 1 & 2 of cycle +# +# #C = y/Min; +# c0 = 0; +# c1 = 1/Min; +# #Ch = x/Max; +# c0h = 0; +# c1h = 1/Max; +# d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); +# D = d1; +# d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); +# Dh = d1h; +# d0 = sum(Qr(i,i+1:j-1)); +# #E = 1-d0-y/Min; +# e0 = 1-d0; +# e1 = -1/Min; +# d0h = sum(Qrh(j,i+1:j-1)); +# #Eh = 1-d0h-x/Max; +# e0h = 1-d0h; +# e1h = -1/Max; +# r1 = Qr(i,i+1)*Qrh(i+1,i); +# R = r1; +# r1h = Qrh(j,j-1)*Qr(j-1,j); +# Rh = r1h; +# +# # For Part 3 of cycle +# +# d3h = sum(Qh(j,i+1:j-1)); +# E3h = 1-d3h; +# d3 = sum(Q(i,i+1:j-1)); +# E3 = 1-d3; +# +# # Define coeficients for equation system +# a0 = -pS+2*pS*Rh-pS*Rh^2+pS*R-2*pS*Rh*R+pS*Rh^2*R; +# a1 = -E3h*Max*c1h*e0*Rh+E3h*Max*c1h*e0; +# a3 = -E3h*Max*c1h*e1*Rh+E3h*Max*c1h*Dh*c1+E3h*Max*c1h*e1+pS*c1h*c1-pS*c1h*c1*Rh; +# +# b0 = -pH+2*pH*R+pH*Rh-2*pH*Rh*R-pH*R^2+pH*Rh*R^2; +# b2 = -Min*E3*e0h*R*c1+Min*E3*e0h*c1; +# b3 = Min*E3*e1h*c1+Min*E3*D*c1h*c1-pH*c1h*c1*R-Min*E3*e1h*R*c1+pH*c1h*c1; +# +# C2 = a3*b2; +# C1 = (-a0*b3+a1*b2+a3*b0); +# C0 = a1*b0; +# # Solve: C2*z^2 + C1*z + C0 = 0 +# z1 = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); +# z2 = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); +# +# # Solution 1 +# x1 = -(b0+b2*z1)/(b3*z1); +# y1 = z1; +# # Solution 2 +# x2 = -(b0+b2*z2)/(b3*z2); +# y2 = z2; +# +# x = x2; +# y = y2; +# +# # fprintf(1,'2nd: i=%d, j=%d: x1=%f, y1=%f, x2=%f, y2=%f\n',i,j,x1,y1,x2,y2); +# +# # Test Standing cycle: assume x=y +# +# C0 = a0; C1 = a1; C2 = a3; +# z1S = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); +# z2S = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); +# +# # Test Hanging cycle: assume x=y +# +# C0 = b0; C1 = b2; C2 = b3; +# z1H = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); +# z2H = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); +# +# # fprintf(1,'2nd: i=%d, j=%d: z1S=%f,: z2S=%f, z1H=%f, z2H=%f\n',i,j,z1S,z2S,z1H,z2H); +# +# else +# +# Eye = EYE(1:j-i-2,1:j-i-2); +# +# # For Part 1 & 2 of cycle +# +# I = i+1:j-2; +# J = i+2:j-1; +# A = Qr(I,J); +# Ah = Qrh(J,I); +# a = Qr(i,J); +# ah = Qrh(j,I); +# b = Qr(I,j); +# bh = Qrh(J,i); +# +# e = 1 - sum(Qr(I,i+2:j),2); +# eh = 1 - sum(Qrh(J,i:j-2),2); +# +# Inv = inv(Eye-A*Ah); +# #C = y/Min + a*Ah*Inv*b; +# c0 = a*Ah*Inv*b; +# c1 = 1/Min; +# #Ch = x/Max + ah*Inv*A*bh; +# c0h = ah*Inv*A*bh; +# c1h = 1/Max; +# d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); +# D = d1+a*eh+a*Ah*Inv*A*eh; +# d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); +# Dh = d1h+ah*Inv*e; +# d0 = sum(Qr(i,i+1:j-1)); +# #E = 1-d0-y/Min+a*Ah*Inv*e; +# e0 = 1-d0+a*Ah*Inv*e; +# e1 = -1/Min; +# d0h = sum(Qrh(j,i+1:j-1)); +# #Eh = 1-d0h-x/Max+ah*Inv*A*eh; +# e0h = 1-d0h+ah*Inv*A*eh; +# e1h = -1/Max; +# r1 = Qr(i,i+1)*Qrh(i+1,i); +# R = r1+a*bh+a*Ah*Inv*A*bh; +# r1h = Qrh(j,j-1)*Qr(j-1,j); +# Rh = r1h+ah*Inv*b; +# +# # For Part 3 of cycle +# +# A3 = Q(I,J); +# A3h = Qh(J,I); +# Inv3 = inv(Eye-A3*A3h); +# +# # For Standing cycle +# d3h = sum(Qh(j,i+1:j-1)); +# c3h = Qh(j,I); +# e3h = 1 - sum(Qh(J,i+1:j-2),2); +# E3h = 1-d3h + c3h*Inv3*A3*e3h; +# +# # For Hanging cycle +# d3 = sum(Q(i,i+1:j-1)); +# c3 = Q(i,J); +# e3 = 1 - sum(Q(I,i+2:j-1),2); +# E3 = 1-d3 + c3*A3h*Inv3*e3; +# +# end +# +# if j-i == 1 # First subdiagonal +# +# if i == 1 +# x = Max; +# y = Max; +# elseif j == n +# x = Min; +# y = Min; +# else +# if pS == 0 +# x = 0; +# y = pH; +# elseif pH == 0 +# x = pS; +# y = 0; +# else +# x = Min*pS/(Min-pH); +# y = Max*pH/(Max-pS); +# end +# end +# +# elseif j-i >= 2 +# if i == 1 +# x = Max*(1-sum(Qh(j,2:j-1))); +# y = Max*(1-sum(Qrh(j,2:j-1))); +# elseif j == n +# x = Min*(1-sum(Q(i,i+1:n-1))); +# y = Min*(1-sum(Qr(i,i+1:n-1))); +# else +# if pS == 0 +# x = 0; +# y = pH; +# elseif pH == 0 +# x = pS; +# y = 0; +# else +# # Define coeficients for equation system +# a0 = pS*c0h*c0+pS*Rh^2*R-2*pS*Rh*R-E3h*Max*c0h*e0*Rh+E3h*Max*c0h*e0+2*pS*Rh+pS*R-pS*c0h*c0*Rh-pS-pS*Rh^2+E3h*Max*c0h*Dh*c0; +# a1 = pS*c1h*c0+E3h*Max*c1h*Dh*c0-E3h*Max*c1h*e0*Rh-pS*c1h*c0*Rh+E3h*Max*c1h*e0; +# a2 = pS*c0h*c1+E3h*Max*c0h*e1-pS*c0h*c1*Rh+E3h*Max*c0h*Dh*c1-E3h*Max*c0h*e1*Rh; +# a3 = -E3h*Max*c1h*e1*Rh+E3h*Max*c1h*Dh*c1+E3h*Max*c1h*e1+pS*c1h*c1-pS*c1h*c1*Rh; +# +# b0 = pH*c0h*c0+pH*Rh*R^2-pH+pH*Rh-2*pH*Rh*R-pH*c0h*c0*R+Min*E3*e0h*c0-Min*E3*e0h*R*c0+Min*E3*D*c0h*c0+2*pH*R-pH*R^2; +# b1 = Min*E3*D*c1h*c0+Min*E3*e1h*c0+pH*c1h*c0-Min*E3*e1h*R*c0-pH*c1h*c0*R; +# b2 = -pH*c0h*c1*R-Min*E3*e0h*R*c1+Min*E3*D*c0h*c1+Min*E3*e0h*c1+pH*c0h*c1; +# b3 = Min*E3*e1h*c1+Min*E3*D*c1h*c1-pH*c1h*c1*R-Min*E3*e1h*R*c1+pH*c1h*c1; +# +# C2 = a2*b3-a3*b2; +# C1 = a0*b3-a1*b2+a2*b1-a3*b0; +# C0 = a0*b1-a1*b0; +# #fprintf(1,'i=%d, j=%d, C0/C2=%f,C1/C2=%f,C2=%f\n',i,j,C0/C2,C1/C2,C2); +# # Solve: C2*z^2 + C1*z + C0 = 0 +# z1 = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); +# z2 = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); +# +# # Solution 1 +# x1 = -(b0+b2*z1)/(b1+b3*z1); +# y1 = z1; +# # Solution 2 +# x2 = -(b0+b2*z2)/(b1+b3*z2); +# y2 = z2; +# +# x = x2; +# y = y2; +# +# # fprintf(1,'End: i=%d, j=%d: x1=%f, y1=%f, x2=%f, y2=%f\n',i,j,x1,y1,x2,y2); +# end +# end +# end +# +# # fprintf(1,'i=%d, j=%d: x=%f, y=%f\n',i,j,x,y); +# +# # min-max +# F(i,j) = x; +# +# # max-min +# F(j,i) = y; +# +# # Fill the transitions matrices +# Q(i,j) = x/Min; +# Qh(j,i) = y/Max; +# Qr(i,j) = y/Min; +# Qrh(j,i) = x/Max; +# +# end +# end +# #fprintf(1,'\n'); +# +# +# T(8,:)=clock; +# +# return F,T + + +if __name__ == "__main__": + from wafo.testing import test_docstrings + test_docstrings(__file__)