# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <codecell>

"""
Created on Thu May 29 11:29:32 2014

@author: chiaracasotto
"""
import numpy as np
import os
import scipy.stats as stat
pi = 3.141592653589793

def simplified_bilinear(T, Gamma, EDPlim, dcroof, EDPvec, RDvec, SPO, bUthd, g, MC):
    
# Estimate fragility parameters for an MDOF system approximated by an equivalent
# SDOF and using the R-mu-T relationships of Ruiz-Garcia & Miranda (2007)

# INPUT
# dy     : roof yield displacement
# EDPlim : median drift values that defines the fragility.
# dcroof : median roof displacement values that defines the fragility.
# buEDP  : dispersion (std of log data) characterizing the lognormal
#          distribution of "limit-state EDP capacity" around EDPlim
# du     : ultimate displacement on the backbone.
# T      : ESDOF period (sec)
# Gamma  : the participation factor for the roof displacement(>1)
# g      : the value of "g" in units compatible with T and dlim, dy.
#          The default is 9.81m/s2, assuming that dy and dlim are in meters.
# N      : Number of realizations to use for Monte Carlo when buthd>0 to
#          incorporate extra uncertainty due to a non-deterministic dlim
#          threshold. If N=0 is provided, a simpler approximation is used
#          instead that is not based on Monte Carlo but on closed form
#          expressions.
# OUTPUT
# Sa50   : the median Sa for the fragility, in units of "g"
# bTSa   : the total dispersion, including capacity and demand dispersions.
# bRSa   : the dispersion due to record-to-record variability (i.e. demand only)
# Sa50cf,bTSacf,bRSacf: Same as above, only provided by the closed-form
#          solution. Only to be used for comparison purposes.

    if T<0.15:
        print 'error: T must be larger than 0.15sec'
        os._exit(1)
    if T>3: print 'T>3 is not advised, but let us proceed'

    dry, du = SPO[0],SPO[-2]
    
    # First translate EDPlim to RD terms and fix non-uniqueness and non-monotonicity.
    # Relationship between EDP and RD defined by pushover analysis
    if len(EDPvec)>2:
        [EDPvec,indy]=np.unique(EDPvec,return_index=True);
        RDvec=RDvec[indy];
    drlim=np.interp(EDPlim,EDPvec,RDvec)

    bUthd = np.array(bUthd)
    if bUthd.any>0 and MC==0:
        print 'error: MC must be different from zero'
        os._exit(1)
    tdrlim = drlim
    # Ductility level mu for each Limit State
    mlim = np.divide(drlim,dry);
    
    # limiting ductility 
    tdrlim[drlim>du] = du  
    tmlim = np.divide(tdrlim,dry) # tmlim is limited by du
    print "mu(LS) = ", mlim

# Simplified, closed-form fragility derivation. We are assuming EDP-RD
    # proportionality (i.e. EDP ~ a*RD + b) so that buthd=buEDP
    c = np.multiply(79.12,np.power(T,1.98))
    R50 = 0.5*(1-c+np.sqrt(np.power(c,2)+c*2*(2*tmlim-1)+1))
     
    # Go to an R value at 85% of the R50 for a biased estimate of "b", the local
    # slope in log-log of the median IDA. Also used for bthd, the record-to-record
    # dispersion.
    Rlim = np.array([np.maximum(0.85*ele,1.00001) for ele in R50])
    Cr50 = 1+np.divide(Rlim-1,c)           
    bthd = 1.957*(1/5.876+1/(11.749*(T+0.1)))*(1-np.exp(-0.73*(Rlim-1)))
    b = 1+np.divide(np.log(Cr50),np.log(Rlim))
    
    # Elastic Limit States
    R50[tmlim<1]= tmlim[tmlim<1]
    b[tmlim<1]=1
    bthd[tmlim<1]=0
    Cr50[tmlim<1]=1
    
    if R50.any>6: print 'warning: R > 6 not advised, proceeding anyway'
    
    Sa50cf = np.divide(np.power(2*pi/T,2),g*Gamma*Cr50)*tdrlim
    bRSacf = bthd/b
    
    # dlim was greater than du. In such cases it makes little sense to
    # include the dispersion of dlim in the total dispersion. The reason is
    # that we are already at collapse, so changes in dlim have no effect at
    # all. The MC approach takes care of this automatically. Obviously this
    # remains a rough binary approximation (include / don't include) for dlim
    # values close to du.
    bTSacf = np.sqrt(np.power(bthd,2) + np.power(bUthd,2))/b       
    bTSacf[drlim>du] = bRSacf[drlim>du]
    
    Say=4*np.power(pi,2)*dry/(g*Gamma*np.power(T,2));
    # find the R50-value corresponding to muf    
    muf = du/dry;
    # Now, use the closed-form estimate to establish bounds on R that are good
    # enough so that you can get the required mu-values given R for the more complex
    # numerical estimation.
    R50f = 0.5*(1-c+np.sqrt(np.power(c,2)+2*c*(2*muf-1)+1))
    b = 1+np.log(1+(R50f-1)/c)/np.log(R50f)
    bthd = 1.957*(1/5.876+1/(11.749*(T+0.1)))*(1-np.exp(-0.739*(R50f-1)))
    bTSaf = np.sqrt(np.power(bthd,2) + np.power(bUthd,2))/b

    # Number of points of R to estimate for interpolating the IDAs
    Nr = 80;
    
    SaT50, bTSa = np.zeros_like(mlim),np.zeros_like(mlim)
    
    if bUthd.all() == float(0):
        # No Montecarlo MC=0, simplified method to account for limit state uncertainty           
    
        # Add a 20% "factor of safety" to make sure you get a large enough maxR
        maxR = [max(np.array([1.5,1.2*min(R50[k]*np.exp(bTSacf[k]),R50f*np.exp(bTSaf[k]))])) for k in range(0,len(R50))]
        R = np.array([np.insert(np.linspace(1,ele,Nr),0,0) for ele in maxR])
        #Define constant relative strength inelastic displacement ratio Cr    
        CR50i = 1 + np.divide(R-1,(79.12*np.power(T,1.98)));
        bCR50i=1.957*(1/5.876+1/(11.749*(T+0.1)))*(1-np.exp(-0.739*(R-1)));
        mu50=R*CR50i;
        mu16=mu50*np.exp(-bCR50i);
        mu84=mu50*np.exp(bCR50i);
        x=[mu16, mu50, mu84];
        
        for i in range(0,len(tmlim)):
            rMC=[np.interp(tmlim[i],ele[i],R[i]) for ele in x]
            SaT50[i] = rMC[1]*Say
            bTSa[i] = (np.log(rMC[0])-np.log(rMC[2]))/2
    else:
   #  Monte Carlo simulation. Note how we use here mlim and not tmlim. The reason
   # is that if mlim is slightly larger tham muf, I still want to have the
   # option of the buEDP uncertainty causing some threshold musample values to
   # become less than muf. So impose the mlim<=mf AFTER you sample to generate
   # musample.
    
        maxR2 = np.array([np.max([1.5,1.2*R50f*np.exp(ele)]) for ele in bTSaf])
        maxR1 = np.array([np.max([1.2,1.2*R50[k]*np.exp(bTSacf[k])]) for k in range(0,len(R50))])
        R = []
        for k in range(0,len(maxR1)):
            if maxR1[k]>maxR2[k]:
                # ultimate point controls
                R.append(np.insert(np.linspace(1,maxR2[k],Nr),0,0))
            else:
                # ultimate point is far enough, provide some denser discretization closer
                # to the median mlim.
                R0 = np.insert(np.linspace(1,maxR1[k],np.round(0.5*Nr)),0,0)
                R.append(np.concatenate((R0,np.linspace(maxR1[k]*1.1,maxR2[k],np.round(0.5*Nr))), axis = 0))
        R = np.array(R)
        
        CR50i = 1 + np.divide(R-1,(79.12*np.power(T,1.98)));
        bCR50i=1.957*(1/5.876+1/(11.749*(T+0.1)))*(1-np.exp(-0.739*(R-1)));
        mu50=R*CR50i;
        mu16=mu50*np.exp(-bCR50i);
        mu84=mu50*np.exp(bCR50i);
        x=[mu16, mu50, mu84];
        
        st = (1./(2.*MC))
        en = (1.-(1./(2.*MC)))
        xp = np.linspace(st,en,MC)
        Sai, Sa = [],[]
        for i in range(0,len(mlim)):
            Sai.append([])
            Sa.append(np.array([]))
            if bUthd[i]>0: #do the sampling otherwise go to the simple approach
                EDPsample = stat.lognorm.ppf(xp,bUthd[i],loc=0,scale=EDPlim[i])
                musample = np.interp(EDPsample,EDPvec,RDvec)/dry;
                #musample = stat.lognorm.ppf(xp,bUthd[i],loc=0,scale=mlim[i])
                musample[musample>muf]=muf
                # Estimate R-values of Sa50 and bRSa that correspond to mlim samples.    
                rMC = [np.interp(musample,ele[i],R[i]) for ele in x]
                allSa50 = [ele*Say for ele in rMC[1]]
                allbSa50 = (np.log(rMC[0])-np.log(rMC[2]))/2
                for j in range(0,MC):
                # for each potential (equiprobable) sample of mlim, get a sample of Sa-values
                # as if we actually had N separate IDA curves, distributed according to the Sa50 and bSa50 values extracted above.
                # In other words,generate N Sa-capacities for N different mlim realizations.
                    if allbSa50[j]>0:
                        realisation = stat.lognorm.ppf(xp,allbSa50[j],loc=0,scale=allSa50[j])
                    else:
                        realisation = np.repeat(allSa50[j],MC)
                    Sai[i].append(realisation)
      
        for i in range(0,len(mlim)):
            # Find median of sampled Sa-values for those drlim with bUthd>0, otherwise got to basic RGM method
            if len(Sai[i])>0:           
                for j in range(1,len(Sai[i])):
                    Sai[i][j] = np.concatenate((Sai[i][j-1],Sai[i][j]))
                Sa[i] = Sai[i][-1]
            
                SaT50[i] = np.median(Sa[i])
                bTSa[i] = np.std(np.log(Sa[i]))
            else:
                # For those bUth that are 0 there is not sampling
                rMC = [np.interp(tmlim[i],ele[i],R[i]) for ele in x]
                SaT50[i] = rMC[1]*Say
                bTSa[i] = (np.log(rMC[0])-np.log(rMC[2]))/2
   
    print "median IM = ", SaT50
    print "total dispersion = ", bTSa
    return [SaT50, bTSa]

