vsldoutlierdetectionrow.c

使用INTEL矢量统计类库的程序,包括以下功能: &#61623 Raw and central moments up to 4th order &#61623 Kurtosis and

/*******************************************************************************
!                             INTEL CONFIDENTIAL
!  Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved.
!  The source code contained  or  described herein and all documents related to
!  the source code ("Material") are owned by Intel Corporation or its suppliers
!  or licensors.  Title to the  Material remains with  Intel Corporation or its
!  suppliers and licensors. The Material contains trade secrets and proprietary
!  and  confidential  information of  Intel or its suppliers and licensors. The
!  Material  is  protected  by  worldwide  copyright  and trade secret laws and
!  treaty  provisions. No part of the Material may be used, copied, reproduced,
!  modified, published, uploaded, posted, transmitted, distributed or disclosed
!  in any way without Intel's prior express written permission.
!  No license  under any  patent, copyright, trade secret or other intellectual
!  property right is granted to or conferred upon you by disclosure or delivery
!  of the Materials,  either expressly, by implication, inducement, estoppel or
!  otherwise.  Any  license  under  such  intellectual property  rights must be
!  express and approved by Intel in writing.
!
!*******************************************************************************
!  Content:
!    Outliers detection  Example Program Text
!******************************************************************************/

#include <stdio.h>

#include "mkl.h"
#include "vsl_ss.h"


/* Uncoment next string for start-up with median's method */
//#define MEDIAN

/* Uncoment next string for using User's weigths*/
#define USER_WEIGHTS


#define DIM 30                   /* dimension of the task */
#define P           DIM	
#define M           3              /* number of BACON algorithm parameters */
#define N           100000         /* number of observations */
#define EPSIILON    2              /* percent of outliers in the observation*/
#define BRNG    VSL_BRNG_MCG31     /* VSL basic generator to be used */
#define SEED    7777777            /* Initial value for stream initialization */

static double TempRandom[N];    /* Vector of uniform variates to form outliers */
static double BaconParams[M];   /* Vector of uniform variates to form outliers */
static double BaconWeights[N];  /* Vector of output BACON coefficients */
static double XX[P][N];         /* Transposed matrix of observations */
static double X[N][P];          /* Input matrix of observations */
static double S[P][P];          /* Covariance matrix of input data */
static double mean[P];          /* Mean vector of input data */
static double meanB[P];         /* Mean vector of "bad" data */

static double OutlierSubset[N][P];
static double OutlierIndexGen[N];

#define RETURN_ON_ERROR                 \
    if(errcode < 0)                     \
    {                                   \
        printf("Error: %i\n", errcode); \
        printf("\nTEST FAILED\n");		\
        return errcode;                 \
    }

main()
{
    MKL_INT n;    /* number of rows(observation) */
    MKL_INT p;    /* number of columns (variables) */
    MKL_INT i, j, k;
    MKL_INT info;
    MKL_INT CounterOutlier, FindOutliers;
    MKL_INT Coincidence,FalseAlarm,LostTarget;
    MKL_INT BaconN;
    MKL_INT CounterWeigth = 0;
    int errcode,errcode_task;    
    MKL_INT storage_format_x;


    double Ro = 1. / ( P * 20 );
    double Epsilon = EPSIILON/100.;

    /* Following variables are used in Cholesky factorization subroutine */
    char uplo;
    double init_method, alpha, beta;

    VSLStreamStatePtr stream_good;
    VSLStreamStatePtr stream_outl;
    
    VSLSSTaskPtr task = 0;

    /* Parameters of the BACON algorithm */
    alpha = 1./N;
    beta  = 0.005;
    BaconN = M;

#ifdef MEDIAN
    init_method = VSL_SS_BACON_MEDIAN_INIT_METHOD;
#endif

#ifndef MEDIAN
    init_method = VSL_SS_BACON_MAHALANOBIS_INIT_METHOD;
#endif
    BaconParams[0] = init_method;
    BaconParams[1] = alpha;
    BaconParams[2] = beta;

    n=N;   p=P;

    storage_format_x   = VSL_SS_MATRIX_ROWS_STORAGE;

    errcode = vsldSSNewTask( &task, &p, &n,  &storage_format_x, (double*)XX,
                                 0, 0 );
    RETURN_ON_ERROR;

    errcode = vsldSSEditOutliersDetection(task, &BaconN,(double *)BaconParams,
        (double *)BaconWeights );
    RETURN_ON_ERROR;

    /***** Printing parametrs *****/
    printf("  Sample of outliers detection.\n");
    printf("--------------------------------------------\n\n");
    printf("Parameters:\n");
    printf("    number of variables  : p = %d\n",p);
    printf("    number of observation: n = %d\n",n);
    printf("                       alpha = %.6f\n",alpha );
    printf("                       beta  = %.6f\n\n",beta);
    
    for( i = 0; i < n; i++ )
    {
        TempRandom[i]   = 0.0;
    }
    
    /* Definition of parameters (covariance matrix, means)
           for input data */
   for ( i = 0; i < p; i++)
   {
       mean[i]  = 0.0;
       meanB[i] = 30.0;

       for( j = 0; j < i; j++ )
       {
           S[i][j] = Ro;
       }

       for( j = i + 1; j < p; j++ )
       {
           S[i][j] = Ro;
       }

       S[i][i] = 1.0;
   }

   /* Here start generation of observations matrix
      (multivarite Gaussian random numbers) */

    uplo = 'U';
    /* MKL Choelsky factorization routine call */
    dpotrf( &uplo, &p, (double*)S, &p, &info );
    if( info != 0 )
    {
        /* Here if function above finished work with error*/
        return  info;
    }

    /* Stream initialization */
    errcode = vslNewStream( &stream_good, BRNG, SEED );
    RETURN_ON_ERROR;

    errcode = vslNewStream( &stream_outl, BRNG, SEED );
    RETURN_ON_ERROR;

    /* Generating random numbers from multivariate normal distribution */
    errcode = vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2,
              stream_good, N, (double *)X, P, VSL_MATRIX_STORAGE_FULL,
              (double*)mean, (double*)S );
    RETURN_ON_ERROR;
  
    /* Generating random numbers from multivariate normal distribution
       to form outliers */
    errcode = vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2,
              stream_outl, N, (double *)OutlierSubset, P,
              VSL_MATRIX_STORAGE_FULL, (double*)meanB, (double*)S );
    RETURN_ON_ERROR;
    
    /* Generating random numbers from uniform distribution to form outliers */
    errcode = vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream_outl,
              N, (double *)TempRandom, 0.0, 1.0 );
    RETURN_ON_ERROR;

    /* Integration of outliers into matrix of observations */
    CounterOutlier = 0;
    for( i = 0; i < n; i++ )
    {
        if( TempRandom[i] < Epsilon )
        {
            CounterOutlier++;
            OutlierIndexGen[i] = 0.0;
            for( k = 0; k < p; k++)
            {
                X[i][k] = OutlierSubset[i][k];
            }
        }
        else
        {
            OutlierIndexGen[i]=1.0;
        }
    }

    for( i = 0; i < n; i++ )
    {
        for( j = 0; j < p; j++ )
        { 
            XX[j][i] = X[i][j];
        }
    }

    /* Stream  deallocation */
    errcode = vslDeleteStream( &stream_good );
    RETURN_ON_ERROR;
    errcode = vslDeleteStream( &stream_outl );
    RETURN_ON_ERROR;
    
    /* Call of BACON function */
    errcode_task = vsldSSCompute( task, VSL_SS_OUTLIERS_DETECTION,
        VSL_SS_BACON_METHOD );

    errcode = vslSSDeleteTask( &task );
    RETURN_ON_ERROR;

    Coincidence  = 0; FindOutliers = 0;
    FalseAlarm   = 0; LostTarget   = 0;

    for( i=0; i < n; i++ )
    {
         if( BaconWeights[i] == 0.0 )
         { 
             FindOutliers++;
             if( OutlierIndexGen[i] == 0.0 )
             {
                 Coincidence++;
             }
             else
             {
                 FalseAlarm++;
             }
         }
         else // BaconWeights[i] != 0.0 
         {
             if( OutlierIndexGen[i] == 0.0 )
             {
                 LostTarget++;
             }
         }
    }
    /***** Printing results *****/
    printf("\n\n\nResults of BACON algorithm\n");
    printf("--------------------------------------------\n\n");
    printf("    number of the generated outliers : %d\n",CounterOutlier);
    printf("     number of the detected outliers : %d\n",FindOutliers);
    printf("\n--------------------------------------------\n\n");
    printf("Check of the detection correctness.\n\n");
    printf("    coincidence = %3.2f%%\n",
          ( CounterOutlier > 0 )? 100.*Coincidence/CounterOutlier : 0);
    printf("    false alarm = %3.2f%%\n",
          ( FindOutliers > 0 )? 100.*FalseAlarm/FindOutliers : 0 );
    printf("    lost target = %3.2f%%\n",
          ( CounterOutlier > 0 )? 100.*LostTarget/CounterOutlier : 0 );
    printf("\n--------------------------------------------\n");

    if(errcode_task != VSL_STATUS_OK )
    {
        printf("\n errcode is  %d\n",errcode_task);
    }
    if( Coincidence == CounterOutlier )
    {
        printf("\nTEST PASSED\n");
    }
    else 
    {
        printf("\nTEST FAILED\n");
    }

    return 0;
}