SUBROUTINE MONTE_CARLO_IMPORTANCE_SIMULATION(BASIC_SAMPLING,RV_SET,JACOB,NLS,DP,SYSTEM_PF,META_MODEL,RESP_SURFACE)

!==============================================================================
!	EVALUATES FAILURE PROBABILITIES BY MONTE CARLO SIMULATION
!==============================================================================
!	Andre Teofilo Beck, Ph.D.						          Created: 02/01/06
!	Rodolfo Krul Tessari        				         Last Updated: 24/08/15
!	Engenharia de Estruturas	EESC / USP
!==============================================================================
    
    USE OPTIMIZATION_VARIABLES_DATA_EXCHANGE
    USE CLASS_MONTE_CARLO_IS
    USE MCS_FILE_MANAGER
    
    USE OPEN_MP_VARIABLES
    
!==============================================================================
!	EXTERNAL VARIABLES:
!
!			   NLS  Number of limit state functions for the problem
!	        RV_SET  Random variable set data
!	         JACOB  Jacobian X to Y transformation object
!				DP	Vector of design point objects
!	     SYSTEM_PF  Statistics of system failure probability
! RESPONSE_SURFACE  Response surface for RS-based simulation
!==============================================================================

	IMPLICIT NONE
	INTEGER									      :: NLS
    TYPE(RANDOM_VARIABLE_SET),      INTENT(IN)    :: RV_SET
    TYPE(JACOBIAN),                 INTENT(INOUT) :: JACOB
    TYPE(DESIGN_POINT),             INTENT(INOUT) :: DP(NLS)
    TYPE(RV_RELIABILITY_RESULTS),   INTENT(INOUT) :: SYSTEM_PF
    TYPE(RESPONSE_SURFACE),OPTIONAL,INTENT(IN)    :: RESP_SURFACE(NLS)
    
    CHARACTER(3),                   INTENT(INOUT) :: BASIC_SAMPLING
    CHARACTER(4),                   INTENT(IN)    :: META_MODEL
    
!==============================================================================
!	INTERNAL VARIABLES:
!
!		     NRV  Number of random variables of the random variable set
!	    NSAMPLES  Internal copy of Number of samples (MC_OPTIONS%NSAMPLES)
!          NMCLS  Number of limit states for MC simulation

!			 NSF  Number of sampling functions
!     REDUCE_NSF  Logical variable, indicates that the number of sampling funcitons should be reduced
!	 REDUCED_NSF  Reduced number of sampling functions

!	   NFAILURES  Cummulative counter of number of failures (NLS+1)
!		FM_ORDER  Order of importance of failure modes
!
!	           G  Vector of limit state function results for each sample (NLS+1)
!	           U  Uniformly distributed random number
!		  PF_SUM  Sum of failure probabilities
!
!	   RV_SAMPLE  Vector of samples, each column is one realization of the set of random variables
!  SAMPLE_WEIGHT  Weight of each individual sample
!      IS_WEIGHT  Weight of each sampling function
!
!	   INDICATOR  Indicator function, contains:
!			  	    ONES where samples correspond to failures and 
!			 	    ZEROS elsewhere
!	  MC_OPTIONS  Options for Monte Carlo simulation
!			  IS  Importance sampling object
!			  SF  Sampling function object
! LIMIT_STATE_PF  Statistics of individual limit state failure probability
!==============================================================================

    INTEGER :: NRV,NSAMPLES,NMCLS
    
    INTEGER :: NSF,REDUCED_NSF
    LOGICAL :: REDUCE_NSF = .TRUE.
    
    INTEGER :: NFAILURES(NLS+1),FM_ORDER(NLS)
    
    REAL(8) :: G(NLS+1),U,PFSUM
    
    REAL(8),ALLOCATABLE :: RV_SAMPLE(:,:),SAMPLE_WEIGHT(:),IS_WEIGHT(:)
    REAL(8),ALLOCATABLE :: SUM1(:),SUM2(:)
    
    INTEGER :: MASTER_INIT
    
    TYPE(CRUDE_MONTE_CARLO_OPTIONS) :: MC_OPTIONS
    TYPE(IMPORTANCE_SAMPLING)    :: IS
    TYPE(SAMPLING_FUNCTION)      :: SF
    TYPE(RV_RELIABILITY_RESULTS) :: LIMIT_STATE_PF
    
    CHARACTER(2) :: PF_ESTIMATION_TECHNIQUE
    
    ! AUXILIARY VARIABLES
    INTEGER,ALLOCATABLE :: NUM_STEP(:)
    INTEGER :: LOCAL_FN_CONV, LOCAL_FN_OUT
    INTEGER :: NSAMPLES_INIT,ANTI_ID=0
    INTEGER :: I,J,K,COUNT,IDEN,NUM,RES
    
    REAL(8) :: VPF(NLS)
    
    LOGICAL :: ALLOC_STATUS
    
    
    !=========================================================================================
    ! 0. INITIALIZATION OF VARIABLES
    !=========================================================================================
    
    PF_ESTIMATION_TECHNIQUE='IS'
    
    NRV = RV_SET%NRV
    NFAILURES = 0
    G = 0.D0; VPF = 0.D0
    
    LSC%COUNTER(1) = 0
    
    IF(NLS.EQ.1) NMCLS = 1
    IF(NLS.GT.1) NMCLS = NLS + 1
    
    ALLOCATE(SUM1(NMCLS),SUM2(NMCLS))
    SUM1 = 0.D0; SUM2 = 0.D0;
    
    !*************************************************************************
    !CALL IS_OPTIONS(ID,LOCAL_FN_CONV,LOCAL_FN_OUT,MC_OPTIONS,ANTI_ID)
    CALL READ_ISRV_SIMULATION_OPTIONS(MC_OPTIONS,FN_IN,"*IMPORTANCE_MC_OPTIONS:")
    
    CALL LOCAL_FILE_IDENTIFICATION(PF_ESTIMATION_TECHNIQUE,NLS,BASIC_SAMPLING,LOCAL_FN_CONV,LOCAL_FN_OUT,ANTI_ID,RESP_SURFACE)
    MC_OPTIONS%BASIC_SAMPLING=BASIC_SAMPLING
    
    CALL ALLOCATE_PF(SYSTEM_PF,LIMIT_STATE_PF,MC_OPTIONS)
    
    IF(MC_OPTIONS%GEN_OPTIONS%CONVERGENCE)THEN
        
        CALL OPEN_CONVERGENCE_FILES(PF_ESTIMATION_TECHNIQUE,NLS,LOCAL_FN_CONV,BASIC_SAMPLING,RESP_SURFACE)
        LIMIT_STATE_PF%NSP=MC_OPTIONS%GEN_OPTIONS%NSP;
             SYSTEM_PF%NSP=MC_OPTIONS%GEN_OPTIONS%NSP;
        WRITE(LOCAL_FN_CONV,11)'N','COV','PF_LOWER','PF','PF_UPPER','BETA_LOWER','BETA','BETA_UPPER'
        
    ELSE
        MC_OPTIONS%GEN_OPTIONS%NSP = 1
                LIMIT_STATE_PF%NSP = 1
                     SYSTEM_PF%NSP = 1
    ENDIF
    
    !=========================================================================================
    ! 1. VERIFIES IF DESIGN POINT INFORMATION IS AVAILABLE - EXITS IF NOT
    !=========================================================================================
    
    DO I=1,NLS
        IF(.NOT.DP(I)%CONVERGENCE) THEN
            CALL PRINT_NONFATAL_ERROR_MESSAGE('DESING POINT NOT AVALIABLE, NUMBER:',I)
            RETURN
        ENDIF
    ENDDO
    
    !=========================================================================================
    ! 2. ASSEMBLES IMPORTANCE SAMPLING VECTORS AND ORDERS FAILURE MODES IN TERMS OF IMPORTANCE 
    !=========================================================================================
    
    NSF = MIN(MC_OPTIONS%NSF,NLS)
    ALLOCATE(IS_WEIGHT(NSF))
    
    PFSUM = 0.D0
    DO I=1,NLS
        PFSUM = PFSUM + DP(I)%FO_PF
    ENDDO
    
    IF(PFSUM .LT. 1.D-100) THEN
        IS_WEIGHT(1:NSF) = 1.D0/DBLE(NSF)
        DO I=1,NLS
            FM_ORDER(I) = I
        ENDDO
    ELSE
        PFSUM = 0.D0
        
        IF(LSC%PARALLEL) THEN !FROM SMALLER TO GREATER PF
            
            CALL FAILURE_MODE_CRESCENT_ORDER(NLS,DP(1:NLS)%FO_PF,FM_ORDER(1:NLS))
            
            DO I=1,NSF; VPF(I) = DP(FM_ORDER(I))%FO_PF; ENDDO
            DO I=1,NSF; PFSUM = PFSUM + 1.D0 - VPF(I);  ENDDO
            
            IF (PFSUM.LT.UUF) THEN; CALL PRINT_NONFATAL_ERROR_MESSAGE('PF SUM NILL, EXITING IMPORTANCE SAMPLING!'); RETURN; ENDIF
            
            DO I=1,NSF; IS_WEIGHT(I) = (1.D0 - VPF(I)) / PFSUM; ENDDO
                
        ELSE !FROM GREATER TO SMALLER PF
            
            CALL FAILURE_MODE_DECRESCENT_ORDER(NLS,DP(1:NLS)%FO_PF,FM_ORDER(1:NLS)) 
				
			DO I=1,NSF; VPF(I)= DP(FM_ORDER(I))%FO_PF; ENDDO
			DO I=1,NSF; PFSUM = PFSUM + VPF(I);        ENDDO
   
			IF (PFSUM.LT.UUF) THEN; CALL PRINT_NONFATAL_ERROR_MESSAGE('PF SUM NILL, EXITING IMPORTANCE SAMPLING!'); RETURN; ENDIF
   
			DO I=1,NSF; IS_WEIGHT(I) = VPF(I) / PFSUM; ENDDO
   
        ENDIF
    ENDIF
	
    REDUCED_NSF = NSF
    
    DO WHILE (IS_WEIGHT(REDUCED_NSF).LT.0.01D0)

		REDUCE_NSF = .TRUE.
		REDUCED_NSF = REDUCED_NSF - 1

		PFSUM = 0.D0

		IF(LSC%PARALLEL) THEN
			DO I=1,REDUCED_NSF; PFSUM = PFSUM + 1 - VPF(I);             ENDDO
			DO I=1,REDUCED_NSF; IS_WEIGHT(I) = (1.D0 - VPF(I)) / PFSUM; ENDDO
		ELSE
			DO I=1,REDUCED_NSF; PFSUM = PFSUM + VPF(I);         ENDDO
			DO I=1,REDUCED_NSF; IS_WEIGHT(I)  = VPF(I) / PFSUM; ENDDO
		ENDIF

		CALL PRINT_WARNING_MESSAGE('FAILURE MODE ELIMINATED FROM SAMPLING FUNCTION: ',I)
        
    ENDDO
    
	NSF = REDUCED_NSF
    
    CALL IMPORTANCE_SAMPLING_ALLOC(NSF,IS)
    CALL SAMPLING_FUNCTION_ALLOC(NRV,NSF,SF)
    
    IS%NSF = NSF
    
    IS%WEIGHT(1:NSF) = IS_WEIGHT(1:NSF)
    
    DO I=1,NSF
        IS%SAMPLED_FM(I) = FM_ORDER(I)
        IS%NSAMPLES(I) = INT(IS%WEIGHT(I)*DBLE(MC_OPTIONS%NSAMPLES))
    ENDDO
    
    NSAMPLES = SUM(IS%NSAMPLES)
    
    !CORRIGINDO O TAMANHO DA AMOSTRA PARA SISTEMAS
    IF(NSAMPLES .NE. MC_OPTIONS%NSAMPLES)THEN
        IS%NSAMPLES(1) = IS%NSAMPLES(1) + ABS(NSAMPLES-MC_OPTIONS%NSAMPLES)
    ENDIF
    
    NSAMPLES = 0
    DO I=1,NSF
        NSAMPLES = NSAMPLES + IS%NSAMPLES(I)
    ENDDO
    
    ALLOCATE (RV_SAMPLE(NRV,NSAMPLES), SAMPLE_WEIGHT(NSAMPLES))
    RV_SAMPLE = 0.D0; SAMPLE_WEIGHT = 0.D0;
    
    
    !=========================================================================================
    ! 3. CONSTRUCTS SAMPLING FUNCTION 
    !=========================================================================================
    
    ! IMPOSES SAME CORRELATION STRUCTURE FOR SAMPLING FUNCTION
    SF%NRV = NRV
    SF%CORRELATED_RV = RV_SET%CORRELATED_RV
    !SF%CORRELATED_RV = .FALSE.
    SF%RHO = RV_SET%RHO
    
    DO I=1,NSF
        
        SF%SAMPLING_SET(1:NRV,I) = RV_SET%RV(1:NRV)
        
        ! IMPOSES DESIGN POINT AS MEAN OF SAMPLING SET:
        SF%SAMPLING_SET(1:NRV,I)%MEAN = DP(FM_ORDER(I))%COORD_X(1:NRV)
        
        DO J=1,NRV
            CALL CRV_PARAMETER_EVAL(SF%SAMPLING_SET(J,I))
            CALL CRV_SET_LIMITS(SF%SAMPLING_SET(J,I))
        ENDDO
    ENDDO
    
    
    !=========================================================================================
    ! 4. CONSTRUCTS SAMPLING FUNCTION AND GENERATES SAMPLES OF THE RANDOM VARIABLE SET
    !=========================================================================================
    U=0.D0
    U = RANDOM_NUMBER(MC_OPTIONS%GEN_OPTIONS%SEED) ! RESETS RANDOM NUMBER GENERATOR
    
    
    !IF USING ACADFRAME
    IF (LSC%LS_STRING(1)(1:13)=='BY_SUBROUTINE')THEN
        IF(FEM_KIND%FEM_STRING(1)=='ACADFRAME')THEN
            CALL STATUS_XV(ALLOC_STATUS)
            IF(.NOT. ALLOC_STATUS) CALL ALLOC_XV(NRV)
            CALL MECHANICAL_SOLVER(0,0)
        ENDIF
    ENDIF
    
    IF(BASIC_SAMPLING=='SIM')THEN
        CALL RV_SET_IMPORTANCE_SAMPLE(NRV,RV_SET,JACOB,IS,SF,NSAMPLES,RV_SAMPLE(1:NRV,1:NSAMPLES),SAMPLE_WEIGHT(1:NSAMPLES),0)
	ELSEIF(BASIC_SAMPLING=='ANT')THEN
        CALL RV_SET_IMPORTANCE_SAMPLE(NRV,RV_SET,JACOB,IS,SF,NSAMPLES,RV_SAMPLE(1:NRV,1:NSAMPLES),SAMPLE_WEIGHT(1:NSAMPLES),1)
    ELSEIF(BASIC_SAMPLING=='LHS')THEN
        CALL RV_SET_ILHS_SAMPLE(NRV,RV_SET,JACOB,IS,SF,NSAMPLES,RV_SAMPLE(1:NRV,1:NSAMPLES),SAMPLE_WEIGHT(1:NSAMPLES),MC_OPTIONS%GEN_OPTIONS%SEED)
    ENDIF
    
    IF(MC_OPTIONS%GEN_OPTIONS%FULL_OUT)THEN
        CALL IS_FULL_FILES(LOCAL_FN_OUT,NLS,BASIC_SAMPLING,RESP_SURFACE)
        DO I=1,NSAMPLES
            WRITE(LOCAL_FN_OUT,"(100G15.5)") ( RV_SAMPLE(K,I), K=1,NRV )
        ENDDO
    ENDIF
    
    NSAMPLES_INIT = MC_OPTIONS%NSAMPLES_INIT
    IF(MC_OPTIONS%GEN_OPTIONS%CONVERGENCE)THEN
        
        ALLOCATE(NUM_STEP(MC_OPTIONS%GEN_OPTIONS%NSP))
        
        NUM=   FLOOR(DBLE(NSAMPLES-NSAMPLES_INIT)/DBLE(MC_OPTIONS%GEN_OPTIONS%NSP-1))
        RES=NINT(MOD(DBLE(NSAMPLES - NSAMPLES_INIT),DBLE(MC_OPTIONS%GEN_OPTIONS%NSP-1)))
        
        NUM_STEP=NUM
        NUM_STEP(1)=NSAMPLES_INIT
        NUM_STEP(MC_OPTIONS%GEN_OPTIONS%NSP)=NUM+RES
        
    ELSE
        
        MC_OPTIONS%GEN_OPTIONS%NSP = 1
        MC_OPTIONS%NSAMPLES_INIT=NSAMPLES
        NSAMPLES_INIT=MC_OPTIONS%NSAMPLES_INIT
        ALLOCATE(NUM_STEP(MC_OPTIONS%GEN_OPTIONS%NSP))
        
        NUM_STEP(1)=NSAMPLES
        
    ENDIF
    
	LIMIT_STATE_PF%NSV(1)=NSAMPLES_INIT
	     SYSTEM_PF%NSV(1)=NSAMPLES_INIT
         
    DO I=2,MC_OPTIONS%GEN_OPTIONS%NSP
        LIMIT_STATE_PF%NSV(I)=LIMIT_STATE_PF%NSV(I-1)+NUM_STEP(I)
             SYSTEM_PF%NSV(I)=LIMIT_STATE_PF%NSV(I)
    ENDDO
    
    
    !=========================================================================================
    ! 5. CALLS LIMIT STATE FUNCTION AND EVALUATES FAILURE FOR EACH FAILURE MODE
    !=========================================================================================
    
    IF(OUTPUT_L(3)) CALL PRINT_PF_STATISTICS_HEAD(FN_OUT,MC_OPTIONS%NSAMPLES,MC_OPTIONS%GEN_OPTIONS%SEED)
    IF(OUTPUT_L(3)) CALL PRINT_IMPORTANCE_SAMPLING_STUFF(FN_OUT,IS)
    
    ! REARRANGE NUM_STEP TO PRINT PARTIAL CONVERGENCE RESULTS
    MASTER_INIT=FLOOR(DBLE(NSAMPLES_INIT/NTHREADS))
    NUM=    FLOOR(DBLE(NSAMPLES/DBLE(NTHREADS) - MASTER_INIT)/DBLE(MC_OPTIONS%GEN_OPTIONS%NSP-1))
    RES=FLOOR(MOD(DBLE(NSAMPLES/DBLE(NTHREADS) - MASTER_INIT),DBLE(MC_OPTIONS%GEN_OPTIONS%NSP-1)))
    
    NUM_STEP=NUM
    NUM_STEP(1)=MASTER_INIT
    NUM_STEP(MC_OPTIONS%GEN_OPTIONS%NSP)=NUM+RES
    
    DO I=2,SIZE(NUM_STEP)
        NUM_STEP(I)=NUM_STEP(I-1)+NUM_STEP(I)
    ENDDO
    
    CALL OMP_SET_NUM_THREADS(1)
    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(IDEN)
    
    COUNT=1
    !$OMP DO SCHEDULE(STATIC)
    DO IDEN=1,NSAMPLES
        
        IF (TID.EQ.0 .AND. IDEN.LE.(NSAMPLES/NTHREADS)) CALL PROGRESS_BAR(DBLE(IDEN)/(NSAMPLES/NTHREADS))
        
        IF(META_MODEL=='NONE') THEN
            
            IF(PRESENT(RESP_SURFACE)) THEN
                G(1:NMCLS) = RESPONSE_SURFACE_FUNCTION(NRV,NLS,NMCLS,RESP_SURFACE(1:NLS),RV_SAMPLE(1:NRV,IDEN),NMCLS)
            ELSE
                G(1:NMCLS) = STRAND_LIMIT_STATE_FUNCTION(NRV,NMCLS,RV_SAMPLE(1:NRV,IDEN),NMCLS)
            ENDIF
            
        ELSEIF(META_MODEL=='NNET') THEN
            
            G(1:NMCLS) = STRAND_NN_LIMIT_STATE_FUNCTION(NRV,RV_SET,JACOB,RV_SAMPLE(1:NRV,IDEN),NN_META_MODEL)
            
        END IF
        
        ! CALCULATE AND UPDATE NUMBER OF FAILURES
        DO I=1,NMCLS
            IF (G(I).LE.0.D0) THEN
                !$OMP CRITICAL (SUM_FAILURES)
                NFAILURES(I) = NFAILURES(I) + 1
                SUM1(I) = SUM1(I) + SAMPLE_WEIGHT(IDEN)
                SUM2(I) = SUM2(I) + SAMPLE_WEIGHT(IDEN)**2
                !$OMP END CRITICAL (SUM_FAILURES)
            ENDIF
        ENDDO
        
        ! EVALUATE SYSTEM FAILURE PROBABILITY AND PRINT PARTIAL CONVERGENCE STATISTICS
        IF(MC_OPTIONS%GEN_OPTIONS%CONVERGENCE .AND. TID.EQ.0)THEN
            IF (MOD(IDEN,NUM_STEP(COUNT)).EQ.0) THEN
                CALL FAST_PFSTATISTICS(COUNT,LOCAL_FN_CONV,SYSTEM_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS,SUM1(NMCLS),SUM2(NMCLS),.TRUE.)
                IF(COUNT.LT.SIZE(NUM_STEP)-1) COUNT=COUNT+1
            ENDIF
        ENDIF
        
        !!!! EVALUATE SYSTEM FAILURE PROBABILITY AND PRINT PARTIAL CONVERGENCE STATISTICS
        !!!DO I=1+COUNT,SIZE(NUM_STEP)
        !!!    IF (MOD(IDEN,NUM_STEP(I)).EQ.0) THEN
        !!!        CALL FAST_PFSTATISTICS(I,LOCAL_FN_CONV,SYSTEM_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS,SUM1(NMCLS),SUM2(NMCLS),.TRUE.)
        !!!        !$OMP ATOMIC
        !!!        COUNT=COUNT+1
        !!!        EXIT
        !!!    ENDIF
        !!!ENDDO
        
    ENDDO !IDEN=1,NSAMPLES
    !$OMP END DO
    
    !$OMP END PARALLEL
    
    
    !=========================================================================================
    ! 6. PRINT SYSTEM FAILURE PROBABILITY STATISTICS
    !=========================================================================================
    
    ! FAILURE PROBABILITY STATISTICS FOR EACH LIMIT STATE FUNCTION
    IF(NLS.GT.1) THEN
        DO I=1,NLS
            
            CALL FAST_PFSTATISTICS(MC_OPTIONS%GEN_OPTIONS%NSP,LOCAL_FN_CONV,LIMIT_STATE_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS,SUM1(I),SUM2(I),.FALSE.)
            
            LIMIT_STATE_PF%PF = LIMIT_STATE_PF%MEAN(LIMIT_STATE_PF%NSP)
            
            CALL CHECK_PROBABILITY('MONTE CARLO STATISTICS ',LIMIT_STATE_PF%PF)
            
            LIMIT_STATE_PF%SD =   LIMIT_STATE_PF%SDEV (LIMIT_STATE_PF%NSP)
            LIMIT_STATE_PF%BETA = LIMIT_STATE_PF%BETAV(LIMIT_STATE_PF%NSP)	
            LIMIT_STATE_PF%LB =   LIMIT_STATE_PF%LBV  (LIMIT_STATE_PF%NSP)
            LIMIT_STATE_PF%UB =   LIMIT_STATE_PF%UBV  (LIMIT_STATE_PF%NSP)
            
            DP(I)%MC_PF = LIMIT_STATE_PF%PF
            
            IF(OUTPUT_L(3)) CALL PRINT_PF_STATISTICS(FN_OUT,I,NFAILURES(I),LIMIT_STATE_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS)
            
        ENDDO
    ENDIF
        
    ! SYSTEM FAILURE PROBABILITY STATISTICS
    CALL FAST_PFSTATISTICS(MC_OPTIONS%GEN_OPTIONS%NSP,LOCAL_FN_CONV,SYSTEM_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS,SUM1(NMCLS),SUM2(NMCLS),.TRUE.)
    
    SYSTEM_PF%PF = SYSTEM_PF%MEAN(SYSTEM_PF%NSP)
    
    CALL CHECK_PROBABILITY('MONTE CARLO STATISTICS ',SYSTEM_PF%PF)
    
    SYSTEM_PF%SD =   SYSTEM_PF%SDEV (SYSTEM_PF%NSP)
    SYSTEM_PF%BETA = SYSTEM_PF%BETAV(SYSTEM_PF%NSP)	
    SYSTEM_PF%LB =   SYSTEM_PF%LBV  (SYSTEM_PF%NSP)
    SYSTEM_PF%UB =   SYSTEM_PF%UBV  (SYSTEM_PF%NSP)
    
    IF(OUTPUT_L(3)) CALL PRINT_PF_STATISTICS(FN_OUT,NMCLS,NFAILURES(NMCLS),SYSTEM_PF,MC_OPTIONS%GEN_OPTIONS%BOUNDS)
    
    IF(NMCLS.EQ.1) DP(1)%MC_PF = SYSTEM_PF%PF
    
    IF(NFAILURES(NMCLS).GE.1) SYSTEM_PF%CONVERGENCE = .TRUE.
    
    ! LAST RESULTS
    IF(PRESENT(RESP_SURFACE)) THEN
        SYSTEM_PF%N_LS_EVAL = SUM(RESP_SURFACE(1:NLS)%NP)
    ELSE
        SYSTEM_PF%N_LS_EVAL = LSC%COUNTER(1)
    ENDIF
    
    CALL PUSH_LIMIT_STATE_COUNTER(LSC%COUNTER(2),LSC%COUNTER(1))
    
    
    !=========================================================================================
    ! 7. FINALIZATION (DEALLOCATION OF DYNAMIC ARRAYS/OBJECTS)
    !=========================================================================================
    CALL PROGRESS_BAR(1.d0)
    
    CALL   SAMPLING_FUNCTION_DEALLOC(SF)
	CALL IMPORTANCE_SAMPLING_DEALLOC(IS)
	CALL DEALLOCATE_PF(SYSTEM_PF,LIMIT_STATE_PF)
    DEALLOCATE (IS_WEIGHT,SAMPLE_WEIGHT,RV_SAMPLE)
    
    
10  FORMAT(A,I10)
11  FORMAT(A20,A20,A20,A20,A20,A20,A20,A20)
    
    RETURN
END SUBROUTINE MONTE_CARLO_IMPORTANCE_SIMULATION