C++ AbstractMatrix类代码示例

//virtual void save(string new_file_name, unsigned long int nvars, unsigned long int nobss, unsigned long int * varindexes, unsigned long int * obsindexes)
SEXP save_R(SEXP New_file_name, SEXP IntPars, SEXP s) {
	CHECK_PTR(s);
	AbstractMatrix * p = (AbstractMatrix*) R_ExternalPtrAddr(s);
	if (p == NULL) {
		error_R("pointer is NULL\n");
		return R_NilValue;
	}

	std::string new_file_name = CHAR(STRING_ELT(New_file_name, 0));
	unsigned long int nvars = (unsigned long int) INTEGER(IntPars)[0];
	unsigned long int nobss = (unsigned long int) INTEGER(IntPars)[1];
	unsigned long int varindexes[nvars];
	unsigned long int obsindexes[nobss];

	for (unsigned long int i = 0; i < nvars; i++)
		varindexes[i] = (unsigned long int) INTEGER(IntPars)[i + 2];
	for (unsigned long int i = 0; i < nobss; i++)
		obsindexes[i] = (unsigned long int) INTEGER(IntPars)[i + 2 + nvars];

	try {
		p->saveAs(new_file_name, nvars, nobss, varindexes, obsindexes);
	} catch (int errcode) {
		error_R("can not save data to file %s\n", new_file_name.c_str());
		return R_NilValue;
	}

	SEXP ret;
	PROTECT(ret = allocVector(LGLSXP, 1));
	LOGICAL(ret)[0] = TRUE;
	UNPROTECT(1);
	return ret;
}

SEXP get_all_varnames_R(SEXP s) {
	CHECK_PTR(s);

	AbstractMatrix * p = (AbstractMatrix*) R_ExternalPtrAddr(s);

	if (p == NULL) {
		error_R("pointer is NULL\n");
		return R_NilValue;
	}

	R_len_t nvars = (R_len_t) 0;

	try {
		nvars = p->getNumVariables();
	} catch (int errcode) {
		return R_NilValue;
	}

	fixedchar tmp;
	SEXP ret;
	PROTECT(ret = allocVector(STRSXP, (R_len_t) nvars));

	try {
		for (unsigned long int i = 0; i < nvars; i++) {
			tmp = p->readVariableName(i);
			SET_STRING_ELT(ret, i, mkChar(tmp.name));
		}
	} catch (int errcode) {
		error_R("something went terribly wrong in get_all_varnames_R\n");
		UNPROTECT(1);
		return ret;
	}
	UNPROTECT(1);
	return ret;
}

SEXP read_variable_float_FileMatrix_R(SEXP nvar, SEXP s) {
	CHECK_PTR(s);
	AbstractMatrix * p = (AbstractMatrix*) R_ExternalPtrAddr(s);
	if (p == NULL) {
		error_R("pointer is NULL\n");
		return R_NilValue;
	}
	unsigned long int nvariable = (unsigned long int) INTEGER(nvar)[0];
	unsigned int nobs = 0;
	try {
		nobs = p->getNumObservations();
	} catch (int errcode) {
		return R_NilValue;
	}
	float * internal_data = new (std::nothrow) float[nobs];

	try {
		p->readVariableAs(nvariable, internal_data);
	} catch (int errcode) {
		return R_NilValue;
	}

	SEXP out;
	PROTECT(out = allocVector(REALSXP, (R_len_t) p->getNumObservations()));
	for (unsigned long int i = 0; i < nobs; i++)
		REAL(out)[i] = (double) internal_data[i];
	UNPROTECT(1);

	delete[] internal_data;

	return out;
}

	SEXP get_nvars_R(SEXP s) {
		//	    cout << "get_nvars_R()" << endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);

		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}

		SEXP out;
		PROTECT(out = allocVector(INTSXP, 1));
		unsigned int nvars = 0;

		try {
			nvars = (unsigned int) p->getNumVariables();
		} catch (int errcode) {
			nvars = 0;
		}

		if (nvars<=0) {
			out = R_NilValue;
		} else {
			INTEGER(out)[0] = nvars;
		}
		UNPROTECT(1);
		return out;
	}

SEXP get_nvars_R(SEXP s) {
	CHECK_PTR(s);

	AbstractMatrix * p = (AbstractMatrix*) R_ExternalPtrAddr(s);

	if (p == NULL) {
		error_R("pointer is NULL\n");
		return R_NilValue;
	}

	SEXP out;
	PROTECT(out = allocVector(INTSXP, 1));
	unsigned int nvars = 0;

	try {
		nvars = (unsigned int) p->getNumVariables();
	} catch (int errcode) {
		nvars = 0;
	}

	if (nvars <= 0) {
		out = R_NilValue;
	} else {
		INTEGER(out)[0] = nvars;
	}
	UNPROTECT(1);
	return out;
}

	SEXP read_variable_double_FileMatrix_R(SEXP nvar, SEXP s) {
		//testDbg << "read_variable_float_FileMatrix_R"<<endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);
		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}
		unsigned long nvariable = (unsigned long) INTEGER(nvar)[0] - 1;
		unsigned int nobs = 0;
		try {
			nobs = p->getNumObservations();
		} catch (int errcode) {
			return R_NilValue;
		}
		double * internal_data = new (std::nothrow) double [nobs];

		try {
			p->readVariableAs(nvariable, internal_data);
		} catch (int errcode) {
			return R_NilValue;
		}

		SEXP out;
		PROTECT(out = allocVector(REALSXP, (R_len_t) p->getNumObservations()));
		for (unsigned long i=0;i< nobs; i++) REAL(out)[i] = internal_data[i];
		delete [] internal_data;

		UNPROTECT(1);

		return out;
	}

	SEXP get_all_obsnames_R(SEXP s) {
		//testDbg << "get_all_obsnames_R"<<endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);

		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}

		//R_len_t nobss = (R_len_t) 0;
		unsigned long int nobss = 0;

		try {
			nobss = p->getNumObservations();
		} catch (int errcode) {
			return R_NilValue;
		}

		FixedChar tmp;
		SEXP ret;
		PROTECT(ret = allocVector(STRSXP, (R_len_t) nobss));

		try {
			for (unsigned long i = 0; i< nobss; i++) {
				tmp = p->readObservationName(i);
				SET_STRING_ELT(ret, i, mkChar(tmp.name));
			}
		} catch (int errcode) {
			error_R("something went terribly wrong in get_all_obsnames_R\n");
			UNPROTECT(1);
			return ret;
		}
		UNPROTECT(1);
		return ret;
	}

	SEXP get_nobs_R(SEXP s) {
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);

		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}

		SEXP out;
		PROTECT(out = allocVector(INTSXP, 1));
		unsigned int nobss = 0;

		try {
			nobss = (unsigned int) p->getNumObservations();
		} catch (int errcode) {
			nobss = 0;
		}

		if (nobss<=0) {
			out = R_NilValue;
		} else {
			INTEGER(out)[0] = nobss;
		}
		UNPROTECT(1);
		return out;
	}

SEXP write_variable_double_FileMatrix_R(SEXP nvar, SEXP data, SEXP s) {
	CHECK_PTR(s);
	AbstractMatrix * p = (AbstractMatrix*) R_ExternalPtrAddr(s);
	if (p == NULL) {
		error_R("pointer is NULL\n");
		return R_NilValue;
	}
	unsigned long int nvariable = (unsigned long int) INTEGER(nvar)[0];
	// here generally should be very careful -- what type of data is IN?

	unsigned int nvars = 0;
	unsigned int nobss = 0;

	try {
		nvars = p->getNumVariables();
	} catch (int errocode) {
		return R_NilValue;
	}

	if (nvariable < 0 || nvariable >= nvars) {
		error_R("nvar (%lu) out of range!\n", nvariable);
		return R_NilValue;
	}

	try {
		nobss = p->getNumObservations();
	} catch (int errcode) {
		return R_NilValue;
	}

	//		float * internal_data = new (std::nothrow) float [nobss];
	double internal_data[nobss];
	if (internal_data == NULL) {
		error_R("internal_data pointer is NULL\n");
		return R_NilValue;
	}

	for (unsigned long int i = 0; i < nobss; i++) {
		internal_data[i] = (double) REAL(data)[i];
	}

	//		Rprintf("\n%lu, %lu\n",nvariable,nobss);
	//		for (unsigned long int i=0;i< nobss;i++) {
	//			Rprintf("%f ",internal_data[i]);
	//		}
	try {
		p->writeVariableAs(nvariable, internal_data);
	} catch (int errcode) {
		error_R("can not write variable %ul\n", nvariable);
	}

	SEXP ret;
	PROTECT(ret = allocVector(LGLSXP, 1));
	LOGICAL(ret)[0] = TRUE;
	UNPROTECT(1);
	return ret;
}

	SEXP write_variable_double_FileMatrix_R(SEXP nvar, SEXP data, SEXP s) {
		//testDbg << "write_variable_double_FileMatrix_R"<<endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);
		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}
		unsigned long nvariable = (unsigned long) INTEGER(nvar)[0] - 1;
		// here generally should be very careful -- what type of data is IN?

		unsigned int nvars = 0;
		unsigned int nobss = 0;

		try {
			nvars = p->getNumVariables();
		} catch (int errocode) {
			return R_NilValue;
		}

		if (nvariable <0 || nvariable >= nvars) {
			error_R("nvar (%lu) out of range!\n",nvariable);
			return R_NilValue;
		}

		try {
			nobss = p->getNumObservations();
		} catch (int errcode) {
			return R_NilValue;
		}

		double * internal_data = new (std::nothrow) double [nobss];

		if (internal_data == NULL) {
			error_R("internal_data pointer is NULL\n");
			return R_NilValue;
		}

		for (unsigned long i=0;i< nobss;i++) {
			internal_data[i] = REAL(data)[i];
		}

		try {
			p->writeVariableAs(nvariable, internal_data);
		} catch (int errcode) {
			delete [] internal_data;
			error_R("can not write variable %ul\n",nvariable);
			return R_NilValue;
		}

		SEXP ret;
		PROTECT(ret = allocVector(LGLSXP, 1));
		LOGICAL(ret)[0] = TRUE;
		delete [] internal_data;

		UNPROTECT(1);
		return ret;
	}

void CorrectnessTest::testRandomReadObservations(){
    testDbg << "testRandomReadObservations" << endl;
    string inputFile = getInputFileName();
    string sumFileName = inputFile + string(".fvf_obssum");
    AbstractMatrix* data = new FileVector ( inputFile, 64 );

    testDbg << "Reading file:" << inputFile << endl;

    unsigned long numVariables = data->getNumVariables();
    unsigned long numObservations = data->getNumObservations();

    double *tmpdat = new (nothrow) double[numVariables];

    testDbg << "Size is " << numVariables << " x " << numObservations << endl;

    int numObservationsToTest = 10;
    int observationIdx[numObservationsToTest];

    unsigned long i;

    TestUtil::initRandomGenerator();
    for (int i=0; i<numObservationsToTest; i++) {
        observationIdx[i] = (rand()*numObservations)/RAND_MAX;
    }

    ifstream sums(sumFileName.c_str());
    testDbg << "Reading sum file: " << sumFileName << endl;

    CPPUNIT_ASSERT(sums.good());

    double *sumData = new double[numVariables];
    for(i=0; i<numObservations; i++) {
        sums >> sumData[i];
    }

    for (i = 0 ; i < numObservationsToTest ; i++ )
	{
	    testDbg << i << "(" << observationIdx[i] << ")" << endl;
	    data->readObservation(observationIdx[i], tmpdat);
	    double calcSumm;

	    calcSumm = summData(tmpdat, numVariables);
	    double relDiff = TestUtil::relativeDifference(calcSumm,sumData[observationIdx[i]]);

	    CPPUNIT_ASSERT(relDiff < 1E-2);
	}

    delete[] tmpdat;
    delete[] sumData;

    delete data;

    testDbg << "Finished" << endl;
}

typename Cholesky<Matrix>::DMatrix Cholesky<Matrix>::doSolve( const AbstractMatrix& b )
{
	if ( this->rowNum() != b.rows() )
	{
		std::cerr<<"The order of matrix is "<<this->rowNum()<<" and the dimension of vector is "<<b.rows()<<std::endl;
		throw COException("Cholesky solving error: the size if not consistent!");
	}
	DMatrix result(b);
	copt_lapack_potrs('U',this->rowNum(),b.cols(),__a,this->lda(),result.dataPtr(),result.lda(),&__info);
	if ( __info != 0 )
		std::cerr<<"Warning in Cholesky solver: solving is wrong!"<<std::endl;
	return result;
}

	//virtual void save(string newFilename, unsigned long nvars, unsigned long nobss, unsigned long * varindexes, unsigned long * obsindexes)
	SEXP save_R(SEXP New_file_name, SEXP IntPars, SEXP s)
	{
		//   dbg<<"save_R"<<endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);
		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}

		string newFilename = CHAR(STRING_ELT(New_file_name,0));
		unsigned long nvars = (unsigned long) INTEGER(IntPars)[0];
		unsigned long nobss = (unsigned long) INTEGER(IntPars)[1];
		unsigned long * varindexes = new (std::nothrow) unsigned long [nvars];
		if (varindexes == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}
		unsigned long * obsindexes = new (std::nothrow) unsigned long [nobss];
		if (obsindexes == NULL) {
			error_R("pointer is NULL\n");
			delete [] varindexes;
			return R_NilValue;
		}

		for (unsigned long i = 0; i < nvars; i++)
			varindexes[i] = (unsigned long) INTEGER(IntPars)[i+2];
		for (unsigned long i = 0; i < nobss; i++) {
			obsindexes[i] = (unsigned long) INTEGER(IntPars)[i+2+nvars];
		}

		try {
			p->saveAs(newFilename,nvars,nobss,varindexes,obsindexes);
		} catch (int errcode) {
			error_R("can not save data to file %s\n",newFilename.c_str());
			delete [] obsindexes;
			delete [] varindexes;
			return R_NilValue;
		}

		SEXP ret;
		PROTECT(ret = allocVector(LGLSXP, 1));
		LOGICAL(ret)[0] = TRUE;
		delete [] obsindexes;
		delete [] varindexes;

		UNPROTECT(1);
		return ret;
	}

	SEXP setReadOnly_R(SEXP s, SEXP readOnly) {
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);

		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}

		bool readonly = LOGICAL(readOnly)[0];

		bool result = p->setReadOnly(readonly);

		SEXP ret;
		PROTECT(ret = allocVector(LGLSXP, 1));
		LOGICAL(ret)[0] = result?TRUE:FALSE;
		UNPROTECT(1);
		return ret;
	}

	// !!!
	SEXP set_cachesizeMb_R(SEXP s, SEXP SizeMB)
	{
		//testDbg << "set_cachesizeMb_R"<<endl;
		AbstractMatrix * p = getAbstractMatrixFromSEXP(s);
		if (p == NULL) {
			error_R("pointer is NULL\n");
			return R_NilValue;
		}
		unsigned long sizeMb = (unsigned long) INTEGER(SizeMB)[0];
		try {
			p->setCacheSizeInMb( sizeMb );
		} catch (int errcode) {
			error_R("cannot reset cache size\n");
			return R_NilValue;
		}

		SEXP ret;PROTECT(ret = allocVector(LGLSXP, 1));LOGICAL(ret)[0] = TRUE;UNPROTECT(1);
		return ret;

	}