C++ linalg::Matrix类代码示例

LinAlg::Matrix<Type> LinAlg::Matrix<Type>::operator ()(unsigned* row_interval, unsigned* column_interval)const
{
    LinAlg::Matrix<Type> Ret;

    if(row_interval[0] < row_interval[1]){
        if(column_interval[0] < column_interval[1]){
            Ret.Init(row_interval[1] - row_interval[0] + 1, column_interval[1] - column_interval[0] + 1);
            for(unsigned i = row_interval[0]; i <= row_interval[1]; ++i)
                for(unsigned j = column_interval[0]; j <= column_interval[1]; ++j)
                    Ret.mat[i - row_interval[0]][j - column_interval[0]] = this->mat[i-1][j-1];

        }else{
            Ret.Init(row_interval[1] - row_interval[0] + 1, column_interval[0] - column_interval[1] + 1);
            for(unsigned i = row_interval[0]; i <= row_interval[1]; ++i)
                for(unsigned j = column_interval[0]; j >= column_interval[1]; --j)
                    Ret.mat[i - row_interval[0]][column_interval[0] - j] = this->mat[i-1][j - 1];
        }

    } else{
        if(column_interval[0] < column_interval[1]){
            Ret.Init(row_interval[0] - row_interval[1] + 1, column_interval[1] - column_interval[0] + 1);
            for(unsigned i = row_interval[0]; i >= row_interval[1]; --i)
                for(unsigned j = column_interval[0]; j <= column_interval[1]; ++j)
                    Ret.mat[row_interval[0] - i][j - column_interval[0]] = this->mat[i-1][j-1];
        }else{
            Ret.Init(row_interval[0] - row_interval[1] + 1, column_interval[0] - column_interval[1] + 1);
            for(unsigned i = row_interval[0]; i >= row_interval[1]; --i)
                for(unsigned j = column_interval[0]; j >= column_interval[1]; --j)
                    Ret.mat[row_interval[0] - i][column_interval[0] - j] = this->mat[i-1][j - 1];
        }
    }

    return Ret;
}

PolynomHandler::Polynom<Type> PolynomHandler::simplify(const Type *num,const Type *den,const unsigned &numSize,const unsigned &denSize)
{
    LinAlg::Matrix<Type> numRoots = Roots<Type>(num, numSize);
    LinAlg::Matrix<Type> denRoots = Roots<Type>(den, denSize);
    unsigned equalRootsCount = 0;

    for(unsigned i = 1; i <= numRoots.getNumberOfRows(); ++i)
          if(PolynomHandler::rootsContainRoot(numRoots(i,1),denRoots))
              ++equalRootsCount;

    LinAlg::Matrix<Type>  newNumRoots(numSize - equalRootsCount,2);
    LinAlg::Matrix<Type>  newDenRoots(denSize - equalRootsCount,2);

    unsigned cont = 1;
    for(unsigned i = 1; i <= numRoots.getNumberOfRows(); ++i)
        if(!PolynomHandler::rootsContainRoot(numRoots(i,1),denRoots))
        {
            newNumRoots(cont,1) = numRoots(i,1);
            newNumRoots(cont,2) = numRoots(i,2);
            cont++;
        }

    cont = 1;
    for(unsigned i = 1; i <= denRoots.getNumberOfRows(); ++i)
        if(!PolynomHandler::rootsContainRoot(denRoots(1,i),numRoots))
        {
            newDenRoots(cont,1) = denRoots(i,1);
            newDenRoots(cont,2) = denRoots(i,2);
            cont++;
        }

    return Polynom<Type>(Root2Poly(newNumRoots),Root2Poly(newDenRoots));

//   return ret;
}

void PolynomHandler::Polynom<Type>::setNum(LinAlg::Matrix<Type> Num)
{
    this->num = initPointer<Type>(Num.getNumberOfColumns());
    this->sizeNum = Num.getNumberOfColumns();

    for (unsigned i = 0; i < Num.getNumberOfColumns(); ++i)
        this->num[i] = Num(1,i+1);
}

void PolynomHandler::Polynom<Type>::setDen(LinAlg::Matrix<Type> Den)
{
    this->den = initPointer<Type>(Den.getNumberOfColumns());
    this->sizeDen = Den.getNumberOfColumns();

    for (unsigned i = 0; i < Den.getNumberOfColumns(); ++i)
        this->den[i] = Den(1,i+1);
}

void OptimizationHandler::RecursiveLeastSquare<Type>::Optimize(LinAlg::Matrix<Type> Input, LinAlg::Matrix<Type> Output)
{
    this->model->setLinearVector(Input, Output(from(1)-->Output.getNumberOfRows(),1));

    LinAlg::Matrix<Type> phi = this->model->getLinearVectorA();
    E = Output(from(1)-->Output.getNumberOfRows(),2) - phi*this->model->getModelCoef();
    K = (P*~phi)/(((phi*P)*~phi) + lambda);
    this->model->setModelCoef(this->model->getModelCoef() + K*E);
    P = (P - (K*(phi*P)))/lambda;
}

LinAlg::Matrix<Type> restrictedOptimizationHandler::activeSet<Type>::activeRestrictions(const LinAlg::Matrix<Type> &A,
                                                                                        const LinAlg::Matrix<Type> &b,
                                                                                        Type tol)
{
    LinAlg::Matrix<Type> restrictionsTest = (A*this->x - b), ind;
    for(unsigned i = 1; i <= restrictionsTest.getNumberOfRows(); ++i)
        if(restrictionsTest(i,1) >= tol)
            ind = ind | Type(i);
    return ind;
}

LinAlg::Matrix<Type> LinAlg::operator~ (LinAlg::Matrix<Type> mat)
{
    LinAlg::Matrix<Type> temp(mat.getNumberOfColumns(), mat.getNumberOfRows());

    for(unsigned i = 1; i <= mat.getNumberOfRows(); i++)
        for(unsigned j = 1; j <= mat.getNumberOfColumns(); j++)
            temp(j, i) = mat(i, j);

    return temp;
}

LinAlg::Matrix<Type> ModelHandler::NFIR<Type>::sim(LinAlg::Matrix<Type> Input)
{
    this->Input  = Input;
    LinAlg::Matrix<Type> TempOutput;

    for(unsigned i = 1; i <= Input.getNumberOfColumns(); ++i){
        this->setLinearVector(Input.GetColumn(i),TempOutput);
        TempOutput = TempOutput|~(this->LinearVectorA*this->ModelCoef);
    }
    this->Output = TempOutput;
    return this->Output;
}

void PolynomHandler::Polynom<Type>::init(LinAlg::Matrix<Type> Num)
{
    this->sizeNum = Num.getNumberOfColumns();
    this->num = initPointer(Num.getNumberOfColumns());
    for (unsigned i = 0; i < Num.getNumberOfColumns(); ++i)
        this->num[i] = (Type) Num(1, i+1);

    this->sizeDen = 1;
    this->den = initPointer<Type>(1);
    this->den[0] = 1;
    this->x = 'x';

}

LinAlg::Matrix<float> SistemasLineares::GaussJacobi(LinAlg::Matrix<float> MatrizUni, unsigned MaxIterations)
{
    //Matriz Resposta
    LinAlg::Matrix<float> MatrizRes(MaxIterations, MatrizUni.getNumberOfColumns());
    LinAlg::Matrix<float> C (MatrizUni.getNumberOfRows(), MatrizUni.getNumberOfColumns() - 1);
    LinAlg::Matrix<float> g (MatrizUni.getNumberOfRows(), 1);
    LinAlg::Matrix<float> x0(C.getNumberOfColumns(), 1);

    //    //Deixa o vetor de chute inicial padronizado como vetor linha
    if(this->X0.getNumberOfColumns() < this->X0.getNumberOfRows())
        ~this->X0;

    //    //Insere o chute inicial na Matriz resposta
    for(unsigned i = 1; i < MatrizRes.getNumberOfColumns() - 1; ++i)
        x0(1,i) = this->X0(1,i);

        //Laço para contar as linhas da MatrizUni e Matriz C.
        for(unsigned i = 1; i <= MatrizUni.getNumberOfRows(); ++i)
        {   //Laço para contar as colunas da MAtrizUni e Matriz C.
            for(unsigned j = 1; j < MatrizUni.getNumberOfColumns(); ++j)
            {
                if(i != j)
                    C(i,j) = - MatrizUni(i,j)/MatrizUni(i,i);//Matriz com a diagonal zerada.
            }
            g(i,1) = MatrizUni(i,MatrizUni.getNumberOfColumns()) / MatrizUni(i,i);//Matriz dos termos independentes.
        }

        MatrizRes = ~x0;
        for(unsigned k = 1; k < MaxIterations; ++k)
        {
            x0 =  (C * x0) + g;
            MatrizRes = MatrizRes || ~x0;
        }
    return MatrizRes;
}

LinAlg::Matrix<float> SistemasLineares::Gauss(LinAlg::Matrix<float> MatrizUni)
{
    LinAlg::Matrix<float> MatrizGauss;
    //Laço para contagem das colunas de MatrizUni.
    for(unsigned i = 1; i < MatrizUni.getNumberOfColumns(); i++)
    {   //Laço para contagem das linhas de MatrizUni.
        for(unsigned j = i + 1; j <= MatrizUni.getNumberOfRows();  j++)
        {
            float m = MatrizUni(j,i)/MatrizUni(i,i);
            //Laço para contagem das colunas da equação.
            for(unsigned z = i ; z <= MatrizUni.getNumberOfColumns(); z++)
                MatrizUni(j,z) = MatrizUni(j,z) - m*MatrizUni(i,z);
        }
    }

    MatrizGauss = LinAlg::Zeros<float>(1, MatrizUni.getNumberOfRows());
    float R;
    for(unsigned i = 1; i <= MatrizUni.getNumberOfRows(); ++i)
    {
        unsigned k = MatrizUni.getNumberOfRows() - i + 1;
        R = 0;
        for(unsigned j = k + 1; j <= MatrizUni.getNumberOfColumns() - 1; ++j)
                R = R + MatrizUni(k, j) * MatrizGauss(1, j);

        MatrizGauss(1, k) = (MatrizUni(k, MatrizUni.getNumberOfColumns()) - R) / MatrizUni(k, k);
    }
    return MatrizGauss;
}

void restrictedOptimizationHandler::activeSet<Type>::KKT(LinAlg::Matrix<Type> A,
                                                         LinAlg::Matrix<Type> b,
                                                         LinAlg::Matrix<Type> &x,
                                                         LinAlg::Matrix<Type> &v)
{
    LinAlg::Matrix<Type> K = (this->QuadMat | (~A)) ||
                             (A | LinAlg::Zeros<Type>(A.getNumberOfRows(),A.getNumberOfRows()));
    LinAlg::Matrix<Type> L = -this->LinMat || b;
    LinAlg::Matrix<Type> S = (((~K)*K)^-1)*(~K)*L;
    x = S(from(1)-->this->LinMat.getNumberOfRows(),1);
    if(A.getNumberOfRows() > 0)
        v = S(from(this->LinMat.getNumberOfRows()+1)-->(this->LinMat.getNumberOfRows() + b.getNumberOfRows()),
              from(1)-->S.getNumberOfColumns());
}

LinAlg::Matrix<Type> PolynomHandler::MultPoly(const LinAlg::Matrix<Type> &lhs, const LinAlg::Matrix<Type> &rhs)
{
    unsigned lhsSize = lhs.getNumberOfColumns();
    unsigned rhsSize = rhs.getNumberOfColumns();
    LinAlg::Matrix<Type>ret(1,lhsSize+rhsSize+1);

    for(unsigned i = 0; i < lhsSize; ++i)
        for(unsigned j = 0; j < rhsSize; ++j)
        {
            ret(1,i+j+1) = ret(1,i+j+1) +  lhs(1,i+1)*rhs(1,j+1);
        }

    return ret;
}

void PlotHandler::plot<Type>::generalPlot(LinAlg::Matrix<Type> X)
{
    customPlot = new QCustomPlot(properties.centralWidget);
    customPlot->setGeometry(QRect(0, 0, this->properties.centralWidget->geometry().width(), this->properties.centralWidget->geometry().height()));
    //customPlot->setGeometry(QRect(this->properties.windowPosX, this->properties.windowPosY, this->properties.windowSizeX, this->properties.windowSizeY));

    // add title layout element:
    if(this->properties.titleFlag)
    {
        customPlot->plotLayout()->insertRow(0);
        customPlot->plotLayout()->addElement(0, 0, new QCPPlotTitle(this->customPlot, this->properties.title.c_str()));
    }

    this->properties.rows = X.getNumberOfRows();
    this->properties.columns = X.getNumberOfColumns();
//        this->setLegend();

    QPen pen;
    // add graphs with different scatter styles:
    for (unsigned i = 0; i < this->properties.rows; ++i)
    {
      customPlot->addGraph();
      pen.setColor(QColor(qSin(i*0.3)*100+100, qSin(i*0.6+0.7)*100+100, qSin(i*0.4+0.6)*100+100));
      // generate data:
      QVector<double> x(this->properties.columns), y(this->properties.columns);
      for (unsigned k = 0; k < this->properties.columns; ++k)
      {
        x[k] = X(i+1,k+1);
        y[k] = k;
      }
      customPlot->graph()->setData(x, y);
      customPlot->graph()->rescaleAxes(true);
      customPlot->graph()->setPen(pen);

      if(this->properties.variablesNameFlag)
          customPlot->graph()->setName("Grafico" + QString::number(i+1));

      customPlot->graph()->setLineStyle(QCPGraph::lsLine);

      if(this->properties.xLabelFlag)
          customPlot->xAxis->setLabel(this->properties.xLabel.c_str());
      if(this->properties.yLabelFlag)
          customPlot->yAxis->setLabel(this->properties.yLabel.c_str());
      // set scatter style:

    }
    customPlot->rescaleAxes();
    customPlot->replot();
    customPlot->repaint();
}

void LinAlg::Balance (LinAlg::Matrix<Type> &matrix_to_balance)
{
    unsigned aux = 0;
    Type radix = FLT_RADIX, sqrdx = radix*radix, s, r, g, f, c;

    while(aux == 0)
    {
        aux = 1;
        for(unsigned i = 1; i <= matrix_to_balance.getNumberOfRows(); i++)
        {
            r = c = 0.0;
            for(unsigned j = 1; j <= matrix_to_balance.getNumberOfColumns(); j++)
                if( j != i)
                {
                    c += std::fabs(matrix_to_balance(j, i));
                    r += std::fabs(matrix_to_balance(i, j));
                }
            if(c && r)
            {
                g = r/radix;
                f = 1.0;
                s = c + r;
                while(c < g)
                {
                    f *= radix;
                    c *= sqrdx;
                }

                g = r*radix;
                while(c > g)
                {
                    f /= radix;
                    c /= sqrdx;
                }
                if((c + r)/f < 0.95*s)
                {
                    aux = 0;
                    g = 1.0/f;
                    for(unsigned j = 1; j <= matrix_to_balance.getNumberOfColumns(); j++)
                    {
                        matrix_to_balance(i, j) *= g;
                        matrix_to_balance(j, i) *= f;
                    }
                }
            }
        }
    }

}