/*
 * skew.c
 *
 * Code generation for function 'skew'
 *
 * C source code generated on: Wed Aug 26 14:59:33 2015
 *
 */

/* Include files */
#include "rt_nonfinite.h"
#include "Select_Ground_Motions.h"
#include "skew.h"
#include "Select_Ground_Motions_emxutil.h"
#include "rdivide.h"
#include "nanmean_.h"
#include "power.h"
#include "isequal.h"
#include "Select_Ground_Motions_mexutil.h"
#include "Select_Ground_Motions_data.h"

/* Variable Definitions */
static emlrtRSInfo im_emlrtRSI = { 44, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRSInfo jm_emlrtRSI = { 45, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRSInfo km_emlrtRSI = { 46, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRSInfo lm_emlrtRSI = { 47, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRTEInfo pc_emlrtRTEI = { 1, 14, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRTEInfo rc_emlrtRTEI = { 44, 1, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRTEInfo sc_emlrtRTEI = { 45, 1, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtRTEInfo tc_emlrtRTEI = { 46, 1, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

static emlrtECInfo mb_emlrtECI = { 2, 44, 6, "skew",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/skew.m"
};

/* Function Definitions */
void eml_error(void)
{
  emlrtPushRtStackR2012b(&jd_emlrtRSI, emlrtRootTLSGlobal);
  emlrtErrorWithMessageIdR2012b(emlrtRootTLSGlobal, &ee_emlrtRTEI,
    "Coder:toolbox:power_domainError", 0);
  emlrtPopRtStackR2012b(&jd_emlrtRSI, emlrtRootTLSGlobal);
}

void skew(const emxArray_real_T *x, emxArray_real_T *s)
{
  int32_T ib;
  int32_T tile[2];
  emxArray_real_T *b_x;
  int32_T i;
  emxArray_real_T *s2;
  real_T n;
  const mxArray *y;
  static const int32_T iv61[2] = { 1, 21 };

  const mxArray *m20;
  char_T cv81[21];
  static const char_T cv82[21] = { 'C', 'o', 'd', 'e', 'r', ':', 'M', 'A', 'T',
    'L', 'A', 'B', ':', 'p', 'm', 'a', 'x', 's', 'i', 'z', 'e' };

  int32_T mv[2];
  emxArray_real_T *x0;
  int32_T iacol;
  boolean_T overflow;
  int32_T jcol;
  boolean_T b15;
  int32_T itilerow;
  emxArray_real_T *b_y;
  emxArray_real_T *m3;
  emxArray_real_T *c_y;
  emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal);

  /* SKEWNESS Skewness. */
  /*    S = SKEWNESS(X) returns the sample skewness of the values in X.  For a */
  /*    vector input, S is the third central moment of X, divided by the cube */
  /*    of its standard deviation.  For a matrix input, S is a row vector */
  /*    containing the sample skewness of each column of X.  For N-D arrays, */
  /*    SKEWNESS operates along the first non-singleton dimension. */
  /*  */
  /*    SKEWNESS(X,0) adjusts the skewness for bias.  SKEWNESS(X,1) is the same */
  /*    as SKEWNESS(X), and does not adjust for bias. */
  /*  */
  /*    SKEWNESS(X,FLAG,DIM) takes the skewness along dimension DIM of X. */
  /*  */
  /*    SKEWNESS treats NaNs as missing values, and removes them. */
  /*  */
  /*    See also MEAN, MOMENT, STD, VAR, KURTOSIS. */
  /*    Copyright 1993-2004 The MathWorks, Inc. */
  /*  The output size for [] is a special case, handle it here. */
  /*  The output size for [] is a special case, handle it here. */
  if (b_isequal(x)) {
    ib = s->size[0] * s->size[1];
    s->size[0] = 1;
    s->size[1] = 1;
    emxEnsureCapacity((emxArray__common *)s, ib, (int32_T)sizeof(real_T),
                      &pc_emlrtRTEI);
    s->data[0] = rtNaN;
  } else {
    /*  Figure out which dimension nanmean will work along. */
    /* dim = find(size(x) ~= 1, 1); */
    /*  Need to tile the output of nanmean to center X. */
    /*  tile = ones(1,2); */
    /*  tile(1) = size(x,1); */
    for (ib = 0; ib < 2; ib++) {
      tile[ib] = 1;
    }

    b_emxInit_real_T(&b_x, 2, &pc_emlrtRTEI, TRUE);
    tile[0] = x->size[0];

    /*  Center X, compute its third and second moments, and compute the */
    /*  uncorrected skewness. */
    emlrtPushRtStackR2012b(&im_emlrtRSI, emlrtRootTLSGlobal);
    ib = b_x->size[0] * b_x->size[1];
    b_x->size[0] = x->size[0];
    b_x->size[1] = x->size[1];
    emxEnsureCapacity((emxArray__common *)b_x, ib, (int32_T)sizeof(real_T),
                      &pc_emlrtRTEI);
    i = x->size[0] * x->size[1];
    for (ib = 0; ib < i; ib++) {
      b_x->data[ib] = x->data[ib];
    }

    b_emxInit_real_T(&s2, 2, &sc_emlrtRTEI, TRUE);
    nanmean_(b_x, s2);
    emlrtPushRtStackR2012b(&sl_emlrtRSI, emlrtRootTLSGlobal);
    emxFree_real_T(&b_x);
    n = 1.0;
    for (i = 0; i < 2; i++) {
      if (tile[i] <= 0) {
        n = 0.0;
      } else {
        n *= (real_T)tile[i];
      }
    }

    if (2.147483647E+9 >= n) {
    } else {
      emlrtPushRtStackR2012b(&uj_emlrtRSI, emlrtRootTLSGlobal);
      y = NULL;
      m20 = mxCreateCharArray(2, iv61);
      for (i = 0; i < 21; i++) {
        cv81[i] = cv82[i];
      }

      emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 21, m20, cv81);
      emlrtAssign(&y, m20);
      error(message(y, &rb_emlrtMCI), &sb_emlrtMCI);
      emlrtPopRtStackR2012b(&uj_emlrtRSI, emlrtRootTLSGlobal);
    }

    emlrtPopRtStackR2012b(&sl_emlrtRSI, emlrtRootTLSGlobal);
    for (ib = 0; ib < 2; ib++) {
      mv[ib] = tile[ib];
    }

    for (ib = 0; ib < 2; ib++) {
      tile[ib] = s2->size[ib] * mv[ib];
    }

    b_emxInit_real_T(&x0, 2, &rc_emlrtRTEI, TRUE);
    ib = x0->size[0] * x0->size[1];
    x0->size[0] = tile[0];
    x0->size[1] = tile[1];
    emxEnsureCapacity((emxArray__common *)x0, ib, (int32_T)sizeof(real_T),
                      &jc_emlrtRTEI);
    if ((tile[0] == 0) || (tile[1] == 0)) {
    } else {
      i = 0;
      ib = 0;
      emlrtPushRtStackR2012b(&tl_emlrtRSI, emlrtRootTLSGlobal);
      emlrtPopRtStackR2012b(&tl_emlrtRSI, emlrtRootTLSGlobal);
      iacol = 0;
      emlrtPushRtStackR2012b(&ul_emlrtRSI, emlrtRootTLSGlobal);
      if (1 > s2->size[1]) {
        overflow = FALSE;
      } else {
        overflow = (s2->size[1] > 2147483646);
      }

      if (overflow) {
        emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
        check_forloop_overflow_error();
        emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
      }

      emlrtPopRtStackR2012b(&ul_emlrtRSI, emlrtRootTLSGlobal);
      for (jcol = 1; jcol <= s2->size[1]; jcol++) {
        emlrtPushRtStackR2012b(&vl_emlrtRSI, emlrtRootTLSGlobal);
        if (1 > mv[0]) {
          b15 = FALSE;
        } else {
          b15 = (mv[0] > 2147483646);
        }

        if (b15) {
          emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
          check_forloop_overflow_error();
          emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
        }

        emlrtPopRtStackR2012b(&vl_emlrtRSI, emlrtRootTLSGlobal);
        for (itilerow = 1; itilerow <= mv[0]; itilerow++) {
          x0->data[ib] = s2->data[iacol];
          i = iacol + 1;
          ib++;
        }

        iacol = i;
      }
    }

    emlrtPopRtStackR2012b(&im_emlrtRSI, emlrtRootTLSGlobal);
    for (ib = 0; ib < 2; ib++) {
      tile[ib] = x->size[ib];
    }

    for (ib = 0; ib < 2; ib++) {
      mv[ib] = x0->size[ib];
    }

    emlrtSizeEqCheck2DFastR2012b(tile, mv, &mb_emlrtECI, emlrtRootTLSGlobal);
    ib = x0->size[0] * x0->size[1];
    x0->size[0] = x->size[0];
    x0->size[1] = x->size[1];
    emxEnsureCapacity((emxArray__common *)x0, ib, (int32_T)sizeof(real_T),
                      &pc_emlrtRTEI);
    i = x->size[0] * x->size[1];
    for (ib = 0; ib < i; ib++) {
      x0->data[ib] = x->data[ib] - x0->data[ib];
    }

    b_emxInit_real_T(&b_y, 2, &pc_emlrtRTEI, TRUE);
    emlrtPushRtStackR2012b(&jm_emlrtRSI, emlrtRootTLSGlobal);
    power(x0, b_y);
    nanmean_(b_y, s2);
    emlrtPopRtStackR2012b(&jm_emlrtRSI, emlrtRootTLSGlobal);

    /*  this is the biased variance estimator */
    emlrtPushRtStackR2012b(&km_emlrtRSI, emlrtRootTLSGlobal);
    for (ib = 0; ib < 2; ib++) {
      tile[ib] = x0->size[ib];
    }

    ib = b_y->size[0] * b_y->size[1];
    b_y->size[0] = tile[0];
    b_y->size[1] = tile[1];
    emxEnsureCapacity((emxArray__common *)b_y, ib, (int32_T)sizeof(real_T),
                      &qc_emlrtRTEI);
    ib = tile[0] * tile[1];
    for (i = 0; i < ib; i++) {
      b_y->data[(int32_T)(1.0 + (real_T)i) - 1] = muDoubleScalarPower(x0->data
        [(int32_T)(1.0 + (real_T)i) - 1], 3.0);
    }

    emxFree_real_T(&x0);
    b_emxInit_real_T(&m3, 2, &tc_emlrtRTEI, TRUE);
    nanmean_(b_y, m3);
    emlrtPopRtStackR2012b(&km_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&lm_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&hd_emlrtRSI, emlrtRootTLSGlobal);
    emxFree_real_T(&b_y);
    for (ib = 0; ib < 2; ib++) {
      tile[ib] = s2->size[ib];
    }

    b_emxInit_real_T(&c_y, 2, &pc_emlrtRTEI, TRUE);
    ib = c_y->size[0] * c_y->size[1];
    c_y->size[0] = 1;
    c_y->size[1] = tile[1];
    emxEnsureCapacity((emxArray__common *)c_y, ib, (int32_T)sizeof(real_T),
                      &qc_emlrtRTEI);
    for (i = 0; i < tile[1]; i++) {
      if (s2->data[i] < 0.0) {
        emlrtPushRtStackR2012b(&id_emlrtRSI, emlrtRootTLSGlobal);
        eml_error();
        emlrtPopRtStackR2012b(&id_emlrtRSI, emlrtRootTLSGlobal);
      }

      c_y->data[i] = muDoubleScalarPower(s2->data[i], 1.5);
    }

    emxFree_real_T(&s2);
    emlrtPopRtStackR2012b(&hd_emlrtRSI, emlrtRootTLSGlobal);
    rdivide(m3, c_y, s);
    emlrtPopRtStackR2012b(&lm_emlrtRSI, emlrtRootTLSGlobal);

    /*  Bias correct the skewness. */
    emxFree_real_T(&c_y);
    emxFree_real_T(&m3);
  }

  emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal);
}

/* End of code generation (skew.c) */