C++ MatrixPtr类代码示例

TEST(Matrix, SparseMatrixTranspose) {
  for (auto height : {10, 50, 100}) {
    for (auto width : {10, 50, 100}) {
      auto nnz = height * width;
      for (auto valueType : {FLOAT_VALUE, NO_VALUE}) {
        for (auto format : {SPARSE_CSR, SPARSE_CSC}) {
          for (auto sparseRate : {0.1, 0.2, 0.5}) {
            MatrixPtr matA = Matrix::createSparseMatrix(
                height, width, size_t(nnz * sparseRate), valueType, format);
            MatrixPtr matB(new CpuSparseMatrix(
                width, height, size_t(nnz * sparseRate), valueType, format));
            matA->randomizeUniform();
            matA->transpose(matB, false);

            /*dense matrix transpose*/
            CpuMatrixPtr matC(new CpuMatrix(height, width));
            matC->copyFrom(*matA);
            MatrixPtr matD(new CpuMatrix(width, height));
            matC->transpose(matD, false);

            /*check result*/
            checkSMatrixEqual2Dense(
                std::dynamic_pointer_cast<CpuSparseMatrix>(matB),
                std::dynamic_pointer_cast<CpuMatrix>(matD));
          }
        }
      }
    }
  }
}

void Reflection::generateReflectionIds()
{
    if (millerCount() == 0)
    {
        std::cout << "Warning! Miller count is 0" << std::endl;
    }
    
    int h = miller(0)->getH();
    int k = miller(0)->getK();
    int l = miller(0)->getL();
    
    cctbx::miller::index<> cctbxMiller = cctbx::miller::index<>(h, k, l);
        for (int i = 0; i < ambiguityCount(); i++)
    {
        MatrixPtr ambiguityMat = matrixForAmbiguity(i);
        cctbx::miller::index<> cctbxTwinnedMiller = ambiguityMat->multiplyIndex(&cctbxMiller);
        
        asym_index asymmetricMiller = asym_index(spaceGroup, asymmetricUnit, cctbxTwinnedMiller);
        
    //    sym_equiv_indices equivMaker = sym_equiv_indices(spaceGroup, cctbxTwinnedMiller);
    //    cctbx::miller::index<> asymmetricMiller = equivMaker(0).h();
       
        int newId = reflectionIdForMiller(asymmetricMiller.h());
    //    int newId = reflectionIdForMiller(cctbxMiller);
        
        reflectionIds.push_back(newId);
    }
}

/** @brief Update by mean vectors and covariance matrices */
TEST_F(TestBCM, CovTest)
{
	// pointer should be NULL initially
	EXPECT_FALSE(m_pSumOfWeightedMeans);
	EXPECT_FALSE(m_pSumOfInvCovs);

	// prediction 1
	update(pMean1, pCov1);
	EXPECT_TRUE(m_pSumOfInvCovs->isApprox(*pSumOfInvCovs1));
	EXPECT_TRUE(m_pSumOfWeightedMeans->isApprox(*pSumOfWeightedMeansByCov1));

	// prediction 2
	update(pMean2, pCov2);
	EXPECT_TRUE(m_pSumOfInvCovs->isApprox(*pSumOfInvCovs2));
	EXPECT_TRUE(m_pSumOfWeightedMeans->isApprox(*pSumOfWeightedMeansByCov2));

	// prediction 3
	update(pMean3, pCov3);
	EXPECT_TRUE(m_pSumOfInvCovs->isApprox(*pSumOfInvCovs3));
	EXPECT_TRUE(m_pSumOfWeightedMeans->isApprox(*pSumOfWeightedMeansByCov3));

	// final
	VectorPtr pMean;
	MatrixPtr pCov;
	get(pMean, pCov);
	EXPECT_TRUE(pCov->isApprox(pCovFinal->diagonal())); // get a variance!!! 
	EXPECT_TRUE(pMean->isApprox(*pMeanByCovFinal));
}

void PowerLayer::forward(PassType passType) {
  Layer::forward(passType);

  MatrixPtr inV0 = getInputValue(0);
  MatrixPtr inV1 = getInputValue(1);

  size_t batchSize = inV1->getHeight();
  size_t dataDim = inV1->getWidth();

  CHECK_EQ(getSize(), dataDim);
  CHECK_EQ(1U, inV0->getWidth());
  CHECK_EQ(batchSize, inV0->getHeight());

  {
    REGISTER_TIMER_INFO("FwResetTimer", getName().c_str());
    reserveOutput(batchSize, dataDim);
  }

  MatrixPtr outV = getOutputValue();

  {
    REGISTER_TIMER_INFO("FwPowerTimer", getName().c_str());
    outV->rowPow(0, *inV1, *inV0);
  }
}

void InterpolationLayer::forward(PassType passType) {
  Layer::forward(passType);

  MatrixPtr weightV = getInputValue(0);
  MatrixPtr inV1 = getInputValue(1);
  MatrixPtr inV2 = getInputValue(2);

  size_t batchSize = inV1->getHeight();
  size_t dataDim = inV1->getWidth();

  CHECK_EQ(dataDim, getSize());
  CHECK_EQ(dataDim, inV2->getWidth());
  CHECK_EQ(batchSize, inV1->getHeight());
  CHECK_EQ(batchSize, inV2->getHeight());

  {
    REGISTER_TIMER_INFO("FwResetTimer", getName().c_str());
    resetOutput(batchSize, dataDim);
  }

  MatrixPtr outV = getOutputValue();

  Matrix::resizeOrCreate(weightLast_, batchSize, 1, false, useGpu_);
  weightLast_->one();
  weightLast_->sub(*weightV);

  REGISTER_TIMER_INFO("FwInterpTimer", getName().c_str());
  // outV = inV1 * weight + inV2 * weightLast
  outV->addRowScale(0, *inV1, *weightV);
  outV->addRowScale(0, *inV2, *weightLast_);
}

void ScaleSubRegionLayer::forward(PassType passType) {
  Layer::forward(passType);
  auto in0 = getInput(0);
  imgH_ = in0.getFrameHeight();
  imgW_ = in0.getFrameWidth();
  if (imgH_ == 0 || imgW_ == 0) {
    auto& conf = config_.inputs(0).scale_sub_region_conf();
    imgH_ = conf.image_conf().img_size_y();
    imgW_ = conf.image_conf().img_size();
  }
  MatrixPtr imgV = in0.value;
  size_t batchSize = imgV->getHeight();
  size_t spatialSize = imgH_ * imgW_;
  channelsNum_ = imgV->getWidth() / spatialSize;
  shape_ = TensorShape({batchSize, channelsNum_, imgH_, imgW_});

  resetOutput(batchSize, imgV->getWidth());
  auto& out = getOutput();
  out.setFrameHeight(imgH_);
  out.setFrameWidth(imgW_);

  MatrixPtr indicesV = getInputValue(1);
  indicesShape_ = TensorShape({batchSize, 6});

  REGISTER_TIMER_INFO("ScaleSubRegionForward", getName().c_str());
  BufferArgs inArgs;
  BufferArgs outArgs;
  inArgs.addArg(*imgV, shape_);
  inArgs.addArg(*indicesV, indicesShape_);
  outArgs.addArg(*out.value, shape_, ASSIGN_TO);
  forward_[0]->calc(inArgs, outArgs);
}

void MultiplexLayer::forward(PassType passType) {
  Layer::forward(passType);

  IVectorPtr copyIds = getInput(0).ids;
  MatrixPtr inV1 = getInputValue(1);
  CHECK_EQ(copyIds->getSize(), inV1->getHeight());
  for (size_t i = 2; i < inputLayers_.size(); i++) {
    CHECK_EQ(inV1->getHeight(), getInputValue(i)->getHeight());
    CHECK_EQ(inV1->getWidth(), getInputValue(i)->getWidth());
  }

  calculateCopySchedule(copyIds, inputLayers_.size() - 1);
  {
    REGISTER_TIMER_INFO("FwResetTimer", getName().c_str());
    reserveOutput(inV1->getHeight(), inV1->getWidth());
  }

  MatrixPtr outV = getOutputValue();
  {
    REGISTER_TIMER_INFO("FwLMultplexingTimer", getName().c_str());
    AsyncGpuBlock block;
    for (const CopyInfo& info : copySchedule_) {
      outV->subMatrix(info.startIdx, info.length, tmpDest_)
          ->copyFrom(*getInputValue(info.copyIdx + 1)
                          ->subMatrix(info.startIdx, info.length, tmpSrc_));
    }
  }

  /* activation */ {
    REGISTER_TIMER_INFO("FwAtvTimer", getName().c_str());
    forwardActivation();
  }
}

void CosSimVecMatLayer::backward(const UpdateCallback& callback) {
  CHECK_EQ(backward_.size(), 1UL) << "Only one forward function needed";

  MatrixPtr inV0 = getInputValue(0);
  MatrixPtr inV1 = getInputValue(1);
  MatrixPtr inG0 = getInputGrad(0);
  MatrixPtr inG1 = getInputGrad(1);
  MatrixPtr outV = getOutputValue();
  MatrixPtr outG = getOutputGrad();

  size_t batchSize = inV0->getHeight();
  CHECK(inV0 && inV1 && inG0 && inG1 && outV && outG);
  REGISTER_TIMER_INFO("BwCosVMTimer", getName().c_str());

  for (size_t i = 0; i < batchSize; i++) {
    tmpRow0->setData(inV0->rowBuf(i));
    tmpRow1->setData(inG0->rowBuf(i));
    tmpMtx0->setData(inV1->rowBuf(i));
    tmpMtx1->setData(inG1->rowBuf(i));
    tmpRow2->setData(outV->rowBuf(i));
    tmpRow3->setData(outG->rowBuf(i));

    BufferArgs inputs;
    BufferArgs outputs;
    inputs.addArg(*tmpRow3);
    inputs.addArg(*tmpRow2);
    inputs.addArg(*tmpMtx0);
    inputs.addArg(*tmpRow0);
    outputs.addArg(*tmpMtx1, ADD_TO);
    outputs.addArg(*tmpRow1, ADD_TO);

    backward_[0]->calc(inputs, outputs);
  }
}

ImagePtr Container::Join() const
{
	const MatrixPtr joined_pixel_data = std::make_shared<Matrix>();

	const unsigned height =
			std::max(
				training_image_->Height(),
				working_image_->Height());

	const unsigned width =
			std::max(
				training_image_->Width(),
				working_image_->Width());

	for (unsigned y = 0; y < height; ++y)
	{
		const PixelsPtr joined_row = std::make_shared<Pixels>();
		for (unsigned x = 0; x < width; ++x)
		{
			if (x < working_image_->Width() && y < working_image_->Height())
			{
				joined_row->push_back(working_image_->PixelData(x, y));
			}
			if (x < training_image_->Width() && y < training_image_->Height())
			{
				joined_row->push_back(training_image_->PixelData(x, y));
			}
		}

		joined_pixel_data->push_back(joined_row);
	}

	return std::make_shared<Image>(joined_pixel_data);
}

/**
 Generates twinned/untwinned reflection IDs for easy searching.
 */
void Reflection::generateReflectionIds()
{
    if (millerCount() == 0)
    {
        std::cout << "Warning! Miller count is 0" << std::endl;
    }
    
    int h = miller(0)->getH();
    int k = miller(0)->getK();
    int l = miller(0)->getL();
    
    vec hkl = new_vector(h, k, l);
    
    for (int i = 0; i < ambiguityCount(); i++)
    {
        MatrixPtr ambiguityMat = matrixForAmbiguity(i);
        
        int asuH, asuK, asuL;
        
        ccp4spg_put_in_asu(spaceGroup, hkl.h, hkl.k, hkl.l, &asuH, &asuK, &asuL);
        vec hklAsu = new_vector(asuH, asuK, asuL);
        
        ambiguityMat->multiplyVector(&hklAsu);
       
        int newId = reflectionIdForMiller(hklAsu);
        
        reflectionIds.push_back(newId);
    }
}

void MKLPackedRecurrentLayer::forwardBatch(int batchSize,
                                           size_t numSequences,
                                           const int* starts) {
  if (!batchValue_) {
    batchValue_.reset(new SequenceToBatch(useGpu_));
  }

  batchValue_->resizeOrCreateBatch(batchSize, numSequences, starts, reversed_);

  batchValue_->copyFromSeq(*output_.value);

  {
    REGISTER_TIMER_INFO("RecurrentFwBatch", getName().c_str());
    /* forward one batch */
    for (size_t n = 0; n < batchValue_->getNumBatch(); n++) {
      MatrixPtr batchValue = batchValue_->getBatchValue(n);

      if (n != 0) {
        MatrixPtr preBatchValue =
            batchValue_->getBatchValue(n - 1, batchValue->getHeight());

        packed_weight_->gemm_compute(preBatchValue, batchValue);
      }
      Argument arg;
      arg.value = batchValue;
      activation_->forward(arg).check();
    }
  }
  batchValue_->copyBackSeq(*output_.value);
}

double Miller::expectedRadius(double spotSize, double mosaicity, vec *hkl)
{
    vec usedHKL;
    
    if (hkl == NULL)
    {
        MatrixPtr newMatrix = MatrixPtr();
        rotateMatrixHKL(latestHRot, latestKRot, 0, matrix, &newMatrix);
        
        usedHKL = new_vector(h, k, l);
        hkl = &usedHKL;
        
        newMatrix->multiplyVector(hkl);
    }
    
    
    spotSize = fabs(spotSize);
    double radMos = fabs(mosaicity) * M_PI / 180;
    
    double distanceFromOrigin = length_of_vector(*hkl);
    
    double spotSizeIncrease = fabs(radMos * distanceFromOrigin);
    
    double radius = (spotSize + spotSizeIncrease);
    
    return radius;
}

void FullGateWriterImpl::write(GatePtr pGate, std::ostream& outputStream) {
	std::string delimeter = "*";
	outputStream << "Gate:";

	unsigned int nbSeq = pGate->getLabelSeq().size();
	for(unsigned int i = 0; i < nbSeq; i++) {
		outputStream  << pGate->getLabelSeq()[i];
		if(i < nbSeq - 1) {
			outputStream << delimeter;
		}
	}

	outputStream << std::endl;
	outputStream << "--Gate cost:" << pGate->getCost();
	outputStream << std::endl;
	outputStream << "--Gate matrix:" << std::endl;

	MatrixPtr pMatrix = pGate->getMatrix();
	int nbRows, nbColumns;
	pMatrix->getSize(nbRows, nbColumns);
	for(int i = 0; i < nbRows; i++) {
		for(int j = 0; j < nbColumns; j++) {
			char printfBuffer[PRINT_BUFFER_LENGTH];
			ComplexVal val = pMatrix->getValue(i,j);
			printVal(printfBuffer, val);
			outputStream << printfBuffer;
		}
		outputStream << std::endl;
	}
}

void ImageTab::calculateAutoThreshold() {
  MatrixPtr matrix = _matrix->selectedMatrix();
  if (matrix) {
    matrix->readLock();
    _lowerThreshold->setText(QString::number(matrix->minValue()));
    _upperThreshold->setText(QString::number(matrix->maxValue()));
    matrix->unlock();
  }
}

void SlopeInterceptLayer::backward(const UpdateCallback& callback) {
  MatrixPtr inG = getInputGrad(0);
  MatrixPtr outG = getOutputGrad();

  if (inG) {
    REGISTER_TIMER_INFO("BwSlopeInterceptTimer", getName().c_str());
    inG->add(*outG, config_.slope());
  }
}