typename PEAddingMethods<Register>::ErrorCode PEAddingMethods<Register>::safe_obfuscate(uint8_t coverage, uint8_t min_len, uint8_t max_len)
{
    PEFile *pe = dynamic_cast<PEFile*>(DAddingMethods<Register>::file);
    if(!pe)
        return ErrorCode::BinaryFileNoPe;

    codePointers.clear();
    relocations.clear();

    QList<uint32_t> fileOffsets;
    ErrorCode ec = getAddressesOffsetsFromTextSection(fileOffsets);
    if(ec != ErrorCode::Success)
        return ec;

    std::uniform_int_distribution<int> prob(0, 99);

    foreach(uint32_t offset, fileOffsets)
    {
        if(prob(DAddingMethods<Register>::r_gen) >= coverage)
            continue;

        BinaryCode<Register> code = generateObfuscationCode(pe->getAddressAtCallInstructionOffset(offset), min_len, max_len);

        uint64_t addr = pe->injectUniqueData(code, codePointers, relocations);
        if(addr == 0)
            return ErrorCode::PeOperationFailed;

        pe->setAddressAtCallInstructionOffset(offset, addr);
    }

    if(!pe->addRelocations(relocations))
        return ErrorCode::PeOperationFailed;

    return ErrorCode::Success;
}

  const char* TaggerImpl::toString() {
    os_.assign("");

#define PRINT \
    for (size_t i = 0; i < x_.size(); ++i) { \
      for (std::vector<const char*>::iterator it = x_[i].begin(); \
           it != x_[i].end(); ++it) \
         os_ << *it << '\t'; \
      os_ << yname(y(i)); \
      if (vlevel_ >= 1) os_ << '/' << prob(i); \
      if (vlevel_ >= 2) { \
        for (size_t j = 0; j < ysize_; ++j)  \
          os_ << '\t' << yname(j) << '/' << prob(i, j); \
      } \
      os_ << '\n'; \
    } \
    os_ << '\n';

    if (nbest_ >= 1) {
      for (size_t n = 0; n < nbest_; ++n) {
        if (!next()) break;
        os_ << "# " << n << " " << prob() << '\n';
        PRINT;
      }
    } else {
      if (vlevel_ >= 1) os_ << "# " << prob() << '\n';
      PRINT;
    }

    return const_cast<const char*>(os_.c_str());

#undef PRINT
  }

 double BM::pdf(uint x,  bool logscale)const{
   if(x>n_)
     return logscale ? BOOM::negative_infinity() : 0;
   if(n_==1){
     double p = x==1 ? prob() : 1-prob();
     return logscale ? log(p) : p;
   }
   return dbinom(x,n_, prob(), logscale);
 }

mat AssetAllocationModel::GetTransitionProbabilities(unsigned int stage) const {
	if (GetStageDependence() == STAGE_INDEPENDENT) {
		return ones(1);
	}
	else if (GetStageDependence() == MARKOV) {
		if (stage == 1) {
			//we suppose root is in non-crisis, fixed state == 1
			mat prob(1, 2);
			prob(0, 0) = 1 - param_.markov_to_crisis_probability;
			prob(0, 1) = param_.markov_to_crisis_probability;
			return prob;
		}
		else {
			//make 2x2 matrix of transitions
			mat prob(2, 2);
			prob(0, 0) = 1 - param_.markov_to_crisis_probability;
			prob(0, 1) = param_.markov_to_crisis_probability;
			prob(1, 0) = param_.markov_from_crisis_probability;
			prob(1, 1) = 1 - param_.markov_from_crisis_probability;
			return prob;
		}
	}
	else {
		throw ScenarioModelException("Unknown dependence model.");
	}
}

void murphy_en(double *me, int deg0, double *dbl_coeff0, int deg1, double *dbl_coeff1, double alpha0, double alpha1, double skewness, int nsm)
{
  int i, j;
  double x, y, sx, sy, theta;
  double v0, v1, xp, e, al0, al1;
  double e0, left, right, theta_left, theta_right, theta_len, theta_inc;
  int interval, nop, nsum;
  int deg[2];
  double *dbl_coeff[2];

  sx=sqrt(assess_area*skewness); sy=sx/skewness;
#ifdef LATTICE
  al0=exp(alpha0)/assess_bound0; /* assuming that special-q in lattice sieving
                                  is of magnitude assess_bound0 */
#else
  al0=exp(alpha0);
#endif
  al1=exp(alpha1);
  deg[0]=deg0; deg[1]=deg1;
  dbl_coeff[0]=dbl_coeff0; dbl_coeff[1]=dbl_coeff1;
  nop=find_optima(deg,dbl_coeff,skewness,&assess_optima);
  qsort(assess_optima,nop,sizeof(double),compare_doubles);

  e=0;
  for (interval=1; interval<nop; interval++) {
    left=assess_optima[interval-1]; right=assess_optima[interval];
    theta_left=atan(left); theta_right=atan(right);
    theta_len=theta_right-theta_left;
    nsum=(int)(((double)nsm)/M_PI*theta_len);
    if (nsum<10) nsum=10;
    theta_inc=theta_len/(double)(nsum);
    e0=0; theta=theta_left+theta_inc/2.;
    for (i=0; i<nsum; i++) {
      theta+=theta_inc;
      x=sx*sin(theta); y=sy*cos(theta);
      v0=dbl_coeff0[0]; xp=1.;
      for (j=1; j<=deg0; j++) { xp*=x; v0*=y; v0+=xp*dbl_coeff0[j]; } 
      v1=dbl_coeff1[0]; xp=1.;
      for (j=1; j<=deg1; j++) { xp*=x; v1*=y; v1+=xp*dbl_coeff1[j]; }
      e0+=prob(al0*v0,assess_bound0)*prob(al1*v1,assess_bound1);
    }
    e0/=nsum;
    e+=(e0*theta_len);
  }
  e/=M_PI;
  *me=e;
  free(assess_optima);
}

Problem read_problem(std::string const path) {
    if(path.empty())
        return Problem(0, 0);
    Problem prob(get_nr_line(path), get_nr_field(path));

    FILE *f = open_c_file(path.c_str(), "r");
    char line[kMaxLineSize];

    uint64_t p = 0;
    for(uint32_t i = 0; fgets(line, kMaxLineSize, f) != nullptr; ++i) {
        char *y_char = strtok(line, " \t");
        float const y = (atoi(y_char)>0)? 1.0f : -1.0f;
        prob.Y[i] = y;
        for(; ; ++p) {
            char *idx_char = strtok(nullptr," \t");
            if(idx_char == nullptr || *idx_char == '\n') break;
            uint32_t idx = static_cast<uint32_t>(atoi(idx_char));
            prob.nr_feature = std::max(prob.nr_feature, idx);
            prob.J[p] = idx-1;
        }
    }

    fclose(f);

    return prob;
}

int Compressor::LoadDecompressionTable(void * comp){
	unsigned short * compressed = (unsigned short *)comp;
	
	// load number of table entries 
	table_entries = ReadShort(compressed)+1;

	// load the values the entries decode to
	for ( int a=0;a<table_entries;a++){
		decoder_table[a].decoded_value = ReadShort(compressed+a+1);
	}
	if ( table_entries == 1){
		return 2*sizeof(unsigned short);
	}

	// store the probability of each entry
	ProbabilityCoder prob(&compressed[table_entries+1],FRACTIONAL_BITS);
	int cp=0;
	for (int a=0;a<table_entries;a++){
		decoder_table[a].cprobability = cp;
		int v = prob.ReadSymbol();
		assert(v);
		if ( v == 0 )
			return -1;
		cp += v;
	}

	assert(cp==(1<<FRACTIONAL_BITS));
	if ( cp != (1<<FRACTIONAL_BITS)){
		return -1;
	}
	decoder_table[table_entries].cprobability = 1<<FRACTIONAL_BITS;
	return (table_entries+1)*sizeof(unsigned short)+prob.GetBytesUsed();
}

void BodyNode::fitWorldTransformParentJointImpl(
    const Eigen::Isometry3d& _target, bool _jointLimit)
{
  Joint* parentJoint = getParentJoint();
  size_t dof = parentJoint->getNumGenCoords();

  if (dof == 0)
    return;

  optimizer::Problem prob(dof);

  // Use the current joint configuration as initial guess
  prob.setInitialGuess(parentJoint->getConfigs());

  // Objective function
  TransformObjFunc obj(this, _target, mSkeleton);
  prob.setObjective(&obj);

  // Joint limit
  if (_jointLimit)
  {
    prob.setLowerBounds(parentJoint->getConfigsMin());
    prob.setUpperBounds(parentJoint->getConfigsMax());
  }

  // Solve with gradient-free local minima algorithm
  optimizer::NloptSolver solver(&prob, NLOPT_LN_BOBYQA);
  solver.solve();

  // Set optimal configuration of the parent joint
  Eigen::VectorXd jointQ = prob.getOptimalSolution();
  parentJoint->setConfigs(jointQ, true, true, true);
}

void brute_iter(context_t *context, task_t *task, struct crypt_data * data_single,
		int (*prob)(struct task_t*, struct context_t*, struct crypt_data*))
{
  int i;
  int count_massive[MAXLENGTH];

  for (i = task->from; i < task->to; i++)
    count_massive[i] = 0;

  for (;;) {
    if (prob (task, context, data_single))
      break;

    for (i = task->to - 1; (i >= task->from) && (count_massive[i] >= context->alph_length -1); --i)
      {
	count_massive[i] = 0;
	task->pass[i] = context->alph[0];
      }

    if (i < task->from)
      break;

    task->pass[i] = context->alph[++count_massive[i]];
  }
}

    void CalibratedModel::calibrate(
        const std::vector<boost::shared_ptr<CalibrationHelper> >& instruments,
        OptimizationMethod& method,
        const EndCriteria& endCriteria,
        const Constraint& additionalConstraint,
        const std::vector<Real>& weights) {

        QL_REQUIRE(weights.empty() ||
                   weights.size() == instruments.size(),
                   "mismatch between number of instruments and weights");

        Constraint c;
        if (additionalConstraint.empty())
            c = *constraint_;
        else
            c = CompositeConstraint(*constraint_,additionalConstraint);
        std::vector<Real> w = weights.empty() ?
                              std::vector<Real>(instruments.size(), 1.0):
                              weights;
        CalibrationFunction f(this, instruments, w);

        Problem prob(f, c, params());
        shortRateEndCriteria_ = method.minimize(prob, endCriteria);
        Array result(prob.currentValue());
        setParams(result);
        Array shortRateProblemValues_ = prob.values(result);

        notifyObservers();
    }

  void LowLevelExplorer::readTransitions
  (istream &stream, unsigned numStates, unsigned numTransitions, RationalParser &parser) {
    unsigned lastFrom(1);
    while (!stream.eof()) {
      string line;
      getline(stream, line);
      if ("" != line) {
	istringstream iss(line);
	unsigned from;
	iss >> from;
	from--;
	unsigned to;
	iss >> to;
	to--;
	string restString;
	getline(iss, restString);
	RationalFunction prob(parser.parseRational(restString));
	if (lastFrom != from) {
	  pmc->finishState();
	  if (lastFrom != from - 1) {
	    throw new runtime_error("Transitions must be ordered by source state.");
	  }
	}
	pmc->addSucc(to, prob);
      }
    }

int main(int argc, char **argv)
{
    if (1 >= argc) {
        fprintf(stderr, "USAGE: mipcl-maxsat [file.cnf|file.wcnf]");
        exit(1);
    }
    char *filename = argv[1];

    try {
        puts("B");
        CMaxSat prob(filename); // 1
	puts("C");
        int nv = read_wcnf(prob, filename);      
	puts("D");
        prob.optimize();
	puts("E");

	if (prob.isSolution()) {
	    double *dpX;
	    int *ipHd;
	    int n = prob.getSolution(dpX, ipHd);
	    // dpX,ipHd	dpX[j] is value of variable whose handle is ipHd[j], j=1,...,n, where n is return value.
            //prob.printSolution("primer.sol"); // (10)
	} else {
	    puts("s UNKNOWN");
	    exit(1);
	}
    }
    catch(CException* pe) {
        std::cerr << pe->getErrorMessage() << std::endl;
        delete pe;
        return 1;
    }
    return 0;
}

// Example useage fisherCosmic2(Photon_Roundness[pho],Photon_Angle[pho])
float NonCollisionBG::fisherCosmic2(float Photon_Roundness, float Photon_Angle){
//======= calculation fisher-discriminant  ===
  TH1F *sigRound          = (TH1F*)file_PHOTON->Get("h_Rho_MC");
  TH1F *sigAngle          = (TH1F*)file_PHOTON->Get("h_Angle_MC");
  TH1F *bkgRound          = (TH1F*)file_COSMIC->Get("h_Rho_cos");
  TH1F *bkgAngle          = (TH1F*)file_COSMIC->Get("h_Angle_cos");

  float sigProb1          = prob(sigRound, Photon_Roundness);
  float bkgProb1          = prob(bkgRound, Photon_Roundness);
  
  float sigProb2          = prob(sigAngle, Photon_Angle);
  float bkgProb2          = prob(bkgAngle, Photon_Angle);
    
  //probability of finding signal or background using shower shape variables(signal like event if prob_2 ~1 and bkg like event if prob_2 ~0)
  return sigProb1*sigProb2/(sigProb1*sigProb2 + bkgProb1*bkgProb2);
}

int main(const int argc, const char **argv)
{
  if (argc == 1) {
    std::shared_ptr<Problem> prob(new SimpleProblem());
    std::unique_ptr<Algorithm> rita(new Ritalgo<SimpleField, SimpleStone>(prob));
    rita->solve();
  }
  else {
    std::unique_ptr<Core> core(new Core(argv[1]));
    const std::string probs = core->get();
    std::shared_ptr<Problem> prob(new ProductionProblem(probs));
    std::unique_ptr<Algorithm> rita(new Ritalgo<SimpleField, SimpleStone>(prob, std::move(core)));
    rita->solve();
  }
  return 0;
}

    void CalibratedModel::calibrate(
                    const vector<shared_ptr<CalibrationHelper> >& instruments,
                    OptimizationMethod& method,
                    const EndCriteria& endCriteria,
                    const Constraint& additionalConstraint,
                    const vector<Real>& weights,
                    const vector<bool>& fixParameters) {

        QL_REQUIRE(weights.empty() || weights.size() == instruments.size(),
                   "mismatch between number of instruments (" <<
                   instruments.size() << ") and weights(" <<
                   weights.size() << ")");

        Constraint c;
        if (additionalConstraint.empty())
            c = *constraint_;
        else
            c = CompositeConstraint(*constraint_,additionalConstraint);
        vector<Real> w =
            weights.empty() ? vector<Real>(instruments.size(), 1.0): weights;

        Array prms = params();
        vector<bool> all(prms.size(), false);
        Projection proj(prms,fixParameters.size()>0 ? fixParameters : all);
        CalibrationFunction f(this,instruments,w,proj);
        ProjectedConstraint pc(c,proj);
        Problem prob(f, pc, proj.project(prms));
        shortRateEndCriteria_ = method.minimize(prob, endCriteria);
        Array result(prob.currentValue());
        setParams(proj.include(result));
        problemValues_ = prob.values(result);

        notifyObservers();
    }

void LVlinear_train(lvError *lvErr, const LVlinear_problem *prob_in, const LVlinear_parameter *param_in, LVlinear_model * model_out){
	try{
		// Input verification: Problem dimensions
		if ((*(prob_in->x))->dimSize != (*(prob_in->y))->dimSize)
			throw LVException(__FILE__, __LINE__, "The problem must have an equal number of labels and feature vectors (x and y).");

		//-- Convert problem
		std::unique_ptr<problem> prob(new problem);
		uint32_t nr_nodes = (*(prob_in->y))->dimSize;
		prob->l = nr_nodes;
		prob->y = (*(prob_in->y))->elt;

		// Create and array of pointers (sparse datastructure)
		std::unique_ptr<feature_node*[]> x(new feature_node*[nr_nodes]);
		prob->x = x.get();

		auto x_in = prob_in->x;
		for (unsigned int i = 0; i < (*x_in)->dimSize; i++){
			// Assign the innermost svm_node array pointers to the array of pointers
			auto xi_in_Hdl = (*x_in)->elt[i];
			x[i] = reinterpret_cast<feature_node*>((*xi_in_Hdl)->elt);
		}

		//-- Convert parameters
		std::unique_ptr<parameter> param(new parameter());
		LVConvertParameter(param_in, param.get());

		// Verify parameters
		const char * param_check = check_parameter(prob.get(), param.get());
		if (param_check != nullptr)
			throw LVException(__FILE__, __LINE__, "Parameter check failed with the following error: " + std::string(param_check));

		// Train model
		model *result = train(prob.get(), param.get());

		// Copy model to LabVIEW memory
		LVConvertModel(result, model_out);

		// Release memory allocated by train
		free_model_content(result);
	}
	catch (LVException &ex) {
		ex.returnError(lvErr);
		// To avoid LabVIEW reading and utilizing bad memory, the dimension sizes of arrays is set to zero
		(*(model_out->label))->dimSize = 0;
		(*(model_out->w))->dimSize = 0;
	}
	catch (std::exception &ex) {
		LVException::returnStdException(lvErr, __FILE__, __LINE__, ex);
		(*(model_out->label))->dimSize = 0;
		(*(model_out->w))->dimSize = 0;
	}
	catch (...) {
		LVException ex(__FILE__, __LINE__, "Unknown exception has occurred");
		ex.returnError(lvErr);
		(*(model_out->label))->dimSize = 0;
		(*(model_out->w))->dimSize = 0;
	}
}

void print_results ()
{
    for (double i = 0; i <= 6; i += .1)
    {
        double curp = prob(i * sqrt (2 / PI));
        printf ("%3.1f\t%.10lf\n", i, curp);
    }
}

/**
 * Write edge by xml writer
 */
void GuiEdge::writeByXmlWriter( xmlTextWriterPtr writer)
{
	EdgeAux::writeByXmlWriter( writer);
	xmlTextWriterWriteAttribute( writer, BAD_CAST "label", BAD_CAST edgeLabel().toAscii().data());
	xmlTextWriterWriteFormatAttribute( writer, BAD_CAST "prob", "%d", prob());
	if (0 != edgeStName().compare ("default", Qt::CaseInsensitive))
		xmlTextWriterWriteAttribute( writer, BAD_CAST "style", BAD_CAST edgeStName().toAscii().data());
}

PRINCIPLE
PLAYER::now ()
{
    for (int i = 0; i < PRINCIPLE_SIZE; ++i)
        if (prob (p[i], PROB_MAX))
            return static_cast < PRINCIPLE > (i);
    return PRINCIPLE_SIZE;
}

 void viterbi_decoding(const vector<unsigned>& sequence, vector<unsigned>& states)
 {
     MatrixXi record_state(sequence.size(), n_state);
     
     MatrixXd prob(1, n_state);
     prob = prior_state;
     
     for(auto pos=0; pos<sequence.size(); ++pos)
     {
         VectorXd prob_tmp(n_state);
         for(auto i=0; i<n_state; ++i)
         {
             double      prob_max = -1;
             unsigned    state_max = 0;
             for(auto j=0; j<n_state; ++j)
             {
                 if (prob(j) * translations(j, i) * obversations(sequence[pos], j) >  prob_max)
                 {
                     prob_max = prob(j) * translations(j, i) * obversations(sequence[pos], j);
                     state_max = j;
                 }
             }
             record_state(pos, i) = state_max;
             prob_tmp(i) = prob_max;
         }
         
         prob = prob_tmp;
     }
     
     unsigned & last = states[sequence.size()-1];
     last = 0;
     for(auto i=0; i<n_state; ++i)
     {
         if (prob(last) < prob(i))
         {
             last = i;
         }
     }
     
     states.resize(sequence.size());
     for(auto pos=sequence.size()-2; pos>=0; ++pos)
     {
         states[pos] = record_state(pos, states[pos+1]);
     }
 }