void getDrivesForArray(const std::string& arrayName,
MDInterface& md,
QueryData& data) {
std::string path(md.getPathByDevName(arrayName));
if (path.empty()) {
LOG(ERROR) << "Could not get file path for " << arrayName;
return;
}
mdu_array_info_t array;
if (!md.getArrayInfo(path, array)) {
return;
}
/* Create a vector of with all expected slot positions. As we work through
* the RAID disks, we remove discovered slots */
std::vector<size_t> missingSlots(array.raid_disks);
std::iota(missingSlots.begin(), missingSlots.end(), 0);
/* Keep track of index in QueryData that have removed slots since we can't
* make safe assumptions about it's original slot position if disk_number >=
* total_disk and we're unable to deteremine total number of missing slots
* until we walk thru all MD_SB_DISKS */
std::vector<size_t> removedSlots;
size_t qdPos = data.size();
for (size_t i = 0; i < MD_SB_DISKS; i++) {
mdu_disk_info_t disk;
disk.number = i;
if (!md.getDiskInfo(path, disk)) {
continue;
}
if (disk.major > 0) {
Row r;
r["md_device_name"] = arrayName;
r["drive_name"] = md.getDevName(disk.major, disk.minor);
r["state"] = getDiskStateStr(disk.state);
if (disk.raid_disk >= 0) {
r["slot"] = INTEGER(disk.raid_disk);
missingSlots.erase(
std::remove(
missingSlots.begin(), missingSlots.end(), disk.raid_disk),
missingSlots.end());
/* We assume that if the disk number is less than the total disk count
* of the array, then it assumes its original slot position; If the
* number is greater than the disk count, then it's not safe to make
* that assumption. We do this check here b/c if a recovery is targeted
* for the same slot, we potentially miss identifying the original slot
* position of the bad disk. */
} else if (disk.raid_disk < 0 && disk.number < array.raid_disks) {
r["slot"] = std::to_string(disk.number);
missingSlots.erase(
std::remove(missingSlots.begin(), missingSlots.end(), disk.number),
missingSlots.end());
/* Mark QueryData position as a removedSlot to handle later*/
} else {
removedSlots.push_back(qdPos);
}
qdPos++;
data.push_back(r);
}
}
/* Handle all missing slots. See `scattered_faulty_and_removed` unit test in
* `./tests/md_tables_tests.cpp`*/
for (const auto& slot : missingSlots) {
if (!removedSlots.empty()) {
data[removedSlots[0]]["slot"] = INTEGER(slot);
removedSlots.erase(removedSlots.begin());
} else {
Row r;
r["md_device_name"] = arrayName;
r["drive_name"] = "unknown";
r["state"] = "removed";
r["slot"] = std::to_string(slot);
data.push_back(r);
}
}
}
SEXP agmart3(SEXP surv2, SEXP score2, SEXP weight2, SEXP strata2,
SEXP sortx, SEXP method2) {
int k, ksave;
int p, istrat, indx2;
double deaths, denom, e_denom;
double hazard, e_hazard, cumhaz;
double temp, time;
double wtsum;
int n, person;
/* pointers to the input data */
double *start, *stop, *event;
double *weight, *score;
int *sort1, *sort2, *strata;
int method; /* integer version of input */
/* output */
SEXP resid2;
double *resid;
n = nrows(surv2);
method = asInteger(method2);
start = REAL(surv2);
stop = start +n;
event = stop +n;
weight= REAL(weight2);
score = REAL(score2);
sort1 = INTEGER(sortx);
sort2 = sort1 + n;
strata= INTEGER(strata2);
PROTECT(resid2 = allocVector(REALSXP, n));
resid = REAL(resid2);
/*
** 'person' walks through the the data from 1 to n,
** sort1[0] points to the largest stop time, sort1[1] the next, ...
** 'time' is a scratch variable holding the time of current interest
** 'indx2' walks through the start times. It will be smaller than
** 'person': if person=27 that means that 27 subjects have stop >=time,
** and are thus potential members of the risk set. If 'indx2' =9,
** that means that 9 subjects have start >=time and thus are NOT part
** of the risk set. (stop > start for each subject guarrantees that
** the 9 are a subset of the 27).
** Basic algorithm: move 'person' forward, adding the new subject into
** the risk set. If this is a new, unique death time, take selected
** old obs out of the sums, add in obs tied at this time, then update
** the cumulative hazard. Everything resets at the end of a stratum.
** The sort order is from large time to small, so we encounter a subject's
** ending time first, then their start time.
** The martingale residual for a subject is
** status - (cumhaz at end of their interval - cumhaz at start)*score
*/
istrat=0;
indx2 =0;
denom =0;
cumhaz =0;
for (person=0; person <n; ) {
p = sort1[person];
if (event[p] ==0) { /* censored */
denom += score[p] * weight[p];
resid[p] = cumhaz * score[p];
person++;
} else {
time = stop[p]; /* found a new, unique death time */
/*
** Remove those subjects whose start time is to the right
** from the risk set, and finish computation of their residual
*/
for (; indx2 <strata[istrat]; indx2++) {
p = sort2[indx2];
if (start[p] < time) break;
denom -= score[p] * weight[p];
resid[p] -= cumhaz * score[p];
}
/*
** Add up over this death time, for all subjects
*/
deaths =0;
e_denom =0;
wtsum =0;
for (k=person; k<strata[istrat]; k++) {
p = sort1[k];
if (stop[p] < time) break; /* only tied times */
denom += score[p] * weight[p];
if (event[p] ==1) {
deaths ++;
e_denom += score[p] * weight[p];
wtsum += weight[p];
}
}
ksave = k;
/* compute the increment in hazard
** hazard = usual increment
** e_hazard = efron increment, for tied deaths only
*/
//.........这里部分代码省略.........
开发者ID:cran，项目名称:survival，代码行数:101，代码来源:agmart3.c
示例4: INTEGER
SEXP gbm_pred
(
SEXP radX, // the data matrix
SEXP rcRows, // number of rows
SEXP rcCols, // number of columns
SEXP rcTrees, // number of trees, may be a vector
SEXP rdInitF, // the initial value
SEXP rTrees, // the list of trees
SEXP rCSplits, // the list of categorical splits
SEXP raiVarType, // indicator of continuous/nominal
SEXP riSingleTree // boolean whether to return only results for one tree
)
{
unsigned long hr = 0;
int iTree = 0;
int iObs = 0;
int cRows = INTEGER(rcRows)[0];
int cPredIterations = LENGTH(rcTrees);
int iPredIteration = 0;
int cTrees = 0;
SEXP rThisTree = NULL;
int *aiSplitVar = NULL;
double *adSplitCode = NULL;
int *aiLeftNode = NULL;
int *aiRightNode = NULL;
int *aiMissingNode = NULL;
int iCurrentNode = 0;
double dX = 0.0;
int iCatSplitIndicator = 0;
bool fSingleTree = (INTEGER(riSingleTree)[0]==1);
SEXP radPredF = NULL;
// allocate the predictions to return
PROTECT(radPredF = allocVector(REALSXP, cRows*cPredIterations));
if(radPredF == NULL)
{
hr = GBM_OUTOFMEMORY;
goto Error;
}
// initialize the predicted values
if(!fSingleTree)
{
// initialize with the intercept for only the smallest rcTrees
for(iObs=0; iObs<cRows; iObs++)
{
REAL(radPredF)[iObs] = REAL(rdInitF)[0];
}
}
else
{
for(iObs=0; iObs<cRows*cPredIterations; iObs++)
{
REAL(radPredF)[iObs] = 0.0;
}
}
iTree = 0;
for(iPredIteration=0; iPredIteration<LENGTH(rcTrees); iPredIteration++)
{
cTrees = INTEGER(rcTrees)[iPredIteration];
if(fSingleTree) iTree=cTrees-1;
if(!fSingleTree && (iPredIteration>0))
{
// copy over from the last rcTrees
for(iObs=0; iObs<cRows; iObs++)
{
REAL(radPredF)[cRows*iPredIteration+iObs] =
REAL(radPredF)[cRows*(iPredIteration-1)+iObs];
}
}
while(iTree<cTrees)
{
rThisTree = VECTOR_ELT(rTrees,iTree);
// these relate to columns returned by pretty.gbm.tree()
aiSplitVar = INTEGER(VECTOR_ELT(rThisTree,0));
adSplitCode = REAL (VECTOR_ELT(rThisTree,1));
aiLeftNode = INTEGER(VECTOR_ELT(rThisTree,2));
aiRightNode = INTEGER(VECTOR_ELT(rThisTree,3));
aiMissingNode = INTEGER(VECTOR_ELT(rThisTree,4));
for(iObs=0; iObs<cRows; iObs++)
{
iCurrentNode = 0;
while(aiSplitVar[iCurrentNode] != -1)
{
dX = REAL(radX)[aiSplitVar[iCurrentNode]*cRows + iObs];
// missing?
if(ISNA(dX))
{
iCurrentNode = aiMissingNode[iCurrentNode];
}
// continuous?
else if(INTEGER(raiVarType)[aiSplitVar[iCurrentNode]] == 0)
{
if(dX < adSplitCode[iCurrentNode])
{
iCurrentNode = aiLeftNode[iCurrentNode];
}
//.........这里部分代码省略.........
SEXP BATCHgdwm_loop_MLPnet (SEXP origNet, SEXP Ptrans, SEXP Ttrans, SEXP nepochs, SEXP rho, SEXP thread_number ) {
//The only difference between wm and without it is the weight update (and the place the values are stored, one is batchgd the other is batchgdwm)
SEXP net;
SEXP R_fcall, args, arg1, arg2, arg3;
PROTECT(net=duplicate(origNet));
int* Ptransdim = INTEGER(coerceVector(getAttrib(Ptrans, R_DimSymbol), INTSXP));
int* Ttransdim = INTEGER(coerceVector(getAttrib(Ttrans, R_DimSymbol), INTSXP));
int n_epochs = INTEGER(nepochs)[0];
struct AMOREnet* ptnet = copynet_RC(net);
struct AMOREneuron** neurons = ptnet->neurons;
/////////////////////////////////////////////////////////////////////////
//Convert input and target to double only once (and instead of copying it every time, just change the pointers)
//Different rows for easy switching pointers
double* input_data = REAL(Ptrans);
double* target_data = REAL(Ttrans);
double** inputs = (double**) R_alloc(Ptransdim[1],sizeof(double*)); //This is an 'Index'
double** targets = (double**) R_alloc(Ptransdim[1],sizeof(double*)); //This is an 'Index'
for (int fila=0; fila < Ptransdim[1]; fila++) {
inputs[fila] = &input_data [fila*Ptransdim[0]];
targets[fila] = &target_data[fila*Ttransdim[0]];
}
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
// Thread number calculation
int n_threads = 1;
#ifdef _OPENMP
{
int max_threads = omp_get_max_threads();
int given_threads = 0;
if (isInteger(thread_number))
given_threads = INTEGER(thread_number)[0];
else if (isNumeric(thread_number))
given_threads = floor(REAL(thread_number)[0]);
if (given_threads <1) //I HAVE THE POWER TO SCHEDULE!
if(max_threads >1)
n_threads = max_threads-1; //Leave a CPU free
else
n_threads = 1;
else if (given_threads > max_threads)
n_threads = max_threads;
else
n_threads = given_threads;
if (neurons[0]->actf == CUSTOM_ACTF) //OMP + R custom functions = bad idea
n_threads = 1;
else if ((ptnet->deltaE.name != LMLS_NAME) && (ptnet->deltaE.name != LMS_NAME))
n_threads = 1;
//printf("Using %i threads from a max of %i.\n",n_threads ,max_threads);
}
#endif
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//Contribution (who is to blame) : Parallelization done by Jose Maria
//Memory allocation for running different threads in parallel:
// Each thread will have his own pool of memory to handle the two kinds of temp vars:
// Vars used only inside the forwards/backwards (v0, v1 and method_delta)
// These vars will be initialized and read only by each thread
// Vars that hold the information on how much the weights and the bias should change
// These vars will be initialized by each thread, then accumulated and read by the master thread when the batch is finished
int n_neurons = ptnet->last_neuron+1;
//Temp values, internal in each iteration
double ** v0s = (double** ) R_alloc(n_threads,sizeof(double* )); //This is an 'Index'
double ** v1s = (double** ) R_alloc(n_threads,sizeof(double* )); //This is an 'Index'
double ** method_deltas = (double** ) R_alloc(n_threads,sizeof(double* )); //This is an 'Index'
//Accumulated values
double ** method_deltas_bias = (double** ) R_alloc(n_threads,sizeof(double* )); //This is an 'Index'
double *** method_sums_delta_x = (double***) R_alloc(n_threads,sizeof(double**)); //This is an 'Index'
for(int id_thread=0; id_thread<n_threads;id_thread++){
double* chunk = (double*) R_alloc(4*n_neurons,sizeof(double)); //Actual chunk of memory for each thread, trying to avoid R_alloc calls
//Advantages: Good proximity reference in cache for values of the same thread, and since it has at least 2 neurons
// (Who would have a NNetwork with less than 2 neurons?), chunks are larger than 64 bytes (i7 L2 cache block size?)
v0s [id_thread] = chunk ;
v1s [id_thread] = &chunk[ n_neurons];
method_deltas [id_thread] = &chunk[2*n_neurons];
method_deltas_bias[id_thread] = &chunk[3*n_neurons];
method_sums_delta_x[id_thread] = (double**) R_alloc(n_neurons,sizeof(double*)); //This is an 'Index'
for(int i=0; i<n_neurons; i++) //Different weigth number for each layer, TODO: R_alloc each layer instead of each neuron
method_sums_delta_x[id_thread][i] = (double*) R_alloc(neurons[i]->last_input_link+1,sizeof(double));
}
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//Consistency (leave pnet as if the function had worked with their values instead of external ones)
// R_alloc should handle freeing the memory, so it's not needed to free the previously allocated memory to avoid memory leaks
// Changing pointer instead of copying data
ptnet->input = inputs[Ptransdim[1]-1];
ptnet->target = targets[Ptransdim[1]-1];
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//.........这里部分代码省略.........
开发者ID:cran，项目名称:AMORE，代码行数:101，代码来源:BATCHgdwm.c
示例13: write_volume
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* Purpose:
*
*
*
*
*
*/
SEXP write_volume(SEXP filename, SEXP nDimensions,
SEXP dimLengths,
SEXP dimStarts,
SEXP dimSteps,
SEXP volumeDataType,
SEXP volumeRange,
SEXP hSlab) {
mihandle_t minc_volume;
midimhandle_t dim[MI2_MAX_VAR_DIMS];
mivolumeprops_t volume_properties;
int result;
int ndx;
int no_dimensions;
int dim_lengths[MI2_MAX_VAR_DIMS];
double dim_starts[MI2_MAX_VAR_DIMS];
double dim_steps[MI2_MAX_VAR_DIMS];
double volume_range_min, volume_range_max;
int volume_data_type;
// pointer to the volume data
double *hSlab_p;
misize_t hSlab_start[MI2_MAX_VAR_DIMS];
misize_t hSlab_count[MI2_MAX_VAR_DIMS];
// start ...
if ( R_DEBUG_mincIO ) Rprintf("write_volume: start ...\n");
/* init number of dimensions and their respective sizes
... yes, I could get this myself, but it's best set in the calling R code */
no_dimensions = INTEGER(nDimensions)[0];
volume_data_type = INTEGER(volumeDataType)[0];
volume_range_min = REAL(volumeRange)[0];
volume_range_max = REAL(volumeRange)[1];
// init volume data pointer
hSlab_p = REAL(hSlab);
// init lenghts, steps, and starts
for (ndx=0; ndx < no_dimensions; ++ndx) {
dim_lengths[ndx] = INTEGER(dimLengths)[ndx];
hSlab_count[ndx] = INTEGER(dimLengths)[ndx];
dim_starts[ndx] = REAL(dimStarts)[ndx];
dim_steps[ndx] = REAL(dimSteps)[ndx];
}
// set the properties for the new output volume
// ... no compression, no chunking, no multi-resolution, ... nothin fancy
result = minew_volume_props(&volume_properties);
if (result != MI_NOERROR) {
error("write_volume:minew_volume_props: Error setting output volume properties: %s.\n", CHAR(STRING_ELT(filename, 0)));
}
result = miset_props_compression_type(volume_properties, MI_COMPRESS_NONE);
if (result != MI_NOERROR) {
error("write_volume:miset_props_compression_type: Error setting output volume properties: %s.\n", CHAR(STRING_ELT(filename, 0)));
}
result = miset_props_multi_resolution(volume_properties, FALSE, 1);
if (result != MI_NOERROR) {
error("write_volume:miset_props_multi_resolution: Error setting output volume properties: %s.\n", CHAR(STRING_ELT(filename, 0)));
}
/* create the 3 output dimensions in the order Z, Y, X, as the volume
is stored in this order */
// z-dim
result = micreate_dimension("zspace",
MI_DIMCLASS_SPATIAL,
MI_DIMATTR_REGULARLY_SAMPLED,
dim_lengths[0],
&dim[0]);
//
if (result != MI_NOERROR) {
error("write_volume: Error initializing the dimension struct for %s dimension.\n", "zspace");
}
// y-dim
result = micreate_dimension("yspace",
MI_DIMCLASS_SPATIAL,
MI_DIMATTR_REGULARLY_SAMPLED,
dim_lengths[1],
&dim[1]);
//
if (result != MI_NOERROR) {
error("write_volume: Error initializing the dimension struct for %s dimension.\n", "yspace");
}
// x-dim
result = micreate_dimension("xspace",
MI_DIMCLASS_SPATIAL,
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论