C++ StateVector::getSize方法代码示例

// -------------------------------------------------------------------------------------------------
//  OperatorList::check(StateVector sv)
//
//  Use: Check the operator list for consitency
//  Input: StateVector sv - the spin configuration
//  Output: bool - true if ok, false if not
// -------------------------------------------------------------------------------------------------
bool OperatorList::check(StateVector sv)
{
    // Check that no diagonal or off-diagonal operator acts
    // on two parallel spins
    int N = sv.getSize();
    int L = getSize();

    for (int i = 0; i < L; i++) {
        int op = list->operators[i].op;
        int spin1 = list->operators[i].spin;
        int spin2 = spin1 + 1;
        if (spin2 > (N - 1))
            spin2 = 0;

        if ((op > 0) && (sv(spin1) == sv(spin2))) {
            cout << "Error: Operating on parallel spins!!" << endl;
            cout << "operator " << op << " on " << spin1 << "," << spin2;
            cout << " in " << sv << endl;
            cout << (*this) << endl;
            return false;
        }

        if (op == 2) {
            sv.flip(spin1);
            sv.flip(spin2);
        }
    }

    return true;
}

/**
 * Copy constructor.
 */
StateVector::StateVector(const StateVector &aVector) :
	_data(0.0)
{
	// INITIAL VALUES
	setNull();

	// SET STATES
	setStates(aVector.getTime(),aVector.getSize(),&aVector.getData()[0]);
}

// -------------------------------------------------------------------------------------------------
//  OperatorList::diagonalMove(StateVector& sv, double beta)
//
//  Use: Performs a diagonal update on the operator list
//  Input: StateVector& sv - the spin configuration
//         double beta - the inverse temperature
//  Output: none
// -------------------------------------------------------------------------------------------------
void OperatorList::diagonalUpdate(StateVector& sv, double beta)
{
    int N = sv.getSize();       // Number of spins
    int L = getSize();          // Operator list length

    for (int i = 0; i < L; i++) {
        int n = countDOD();   // Number of diagonal and off-diagonal operators
        int op = list->operators[i].op;
        if (op == 0) {
            // unit operator
            // Calculate probability to change operator
            double P = (double)(N*beta/(2*(L-n)));
            // Generate a random spin that this operator acts on
            int spin1 = (int)(ran()*N);
            // and its adjacent spin
            int spin2 = (spin1 + 1) % N;
            
            // Change operator to diagonal only if the spins are _not_ parallel...
            if (sv(spin1) != sv(spin2)) {
                // ...and with probability P
                double r = ran();

                if (r < P) {
                    //cout << "Change..." << endl;
                    list->operators[i].op = 1;
                    list->operators[i].spin = spin1;
                }
            }
        } else if (op == 1) { 
            // diagonal operator
            // Calculate probability to change to unit operator
            double P = (double)(2*(L-n+1)/(N*beta));
            double r = ran();
            
            // Get the spins this off-diagonal operator acts on
            int spin1 = list->operators[i].spin;
            int spin2 = (spin1 + 1) % N;
            
            if (r < P) {
                list->operators[i].op = 0;    
                // unit ops cannot act on spins
                list->operators[i].spin = -1;    
            }
        } else if (op == 2) { 
            // off-diagonal operator
            // Get the spins this off-diagonal operator acts on
            int spin1 = list->operators[i].spin;
            int spin2 = (spin1 + 1) % N;
            
            // Flip both spins
            sv.flip(spin1);
            sv.flip(spin2);                
        }
    }
} 

void testGait2354()
{
	ForwardTool forward("subject01_Setup_Forward.xml");
	forward.run();
	Storage results("Results/subject01_states.sto");
	Storage* standard = new Storage();
    string statesFileName("std_subject01_walk1_states.sto");
    forward.loadStatesStorage( statesFileName, standard );

	Array<double> data;
	StateVector* state = results.getStateVector(0);
	double time = state->getTime();
	data.setSize(state->getSize());
	standard->getDataAtTime(time, state->getSize(), data);

	// At initial time all states should be identical except for locked joints may vary slightly due to 
	// differences in OpenSim's integrator and SimTK's
	int nc = forward.getModel().getNumCoordinates();
	for (int j = 0; j < state->getSize(); ++j) {
		stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard->getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j];
        ASSERT_EQUAL(data[j], state->getData()[j], 1e-4, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str() << endl;;
	}

	int i = results.getSize()-1;
	state = results.getStateVector(i);
	time = state->getTime();
	standard->getDataAtTime(time, state->getSize(), data);

	// NOTE: Gait model is running forward open-loop. We cannot expect all the states to
	// be "bang on" and we expect a gradual drift in the coordinates.  Check to see that
	// coordinates haven't drifted too far off.
	for (int j = 0; j < nc; ++j) {
	    stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard->getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j];
        ASSERT_EQUAL(data[j], state->getData()[j], 4e-2, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str() << endl;
	}
}

void testArm26() {
    ForwardTool forward("arm26_Setup_Forward.xml");
	forward.run();

	Storage results("Results/arm26_states.sto");
	Storage* standard = new Storage();
    string statesFileName("std_arm26_states.sto");
    forward.loadStatesStorage( statesFileName, standard );
	StateVector* state = results.getStateVector(0);
	double time = state->getTime();
	Array<double> data;
	data.setSize(state->getSize());
	standard->getDataAtTime(time, state->getSize(), data);
	for (int j = 0; j < state->getSize(); ++j) {
		stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard->getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j] << endl;
        ASSERT_EQUAL(data[j], state->getData()[j], 1.0e-3, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str();
    }

	int i = results.getSize()-1;
    state = results.getStateVector(i);
    time = state->getTime();
	data.setSize(state->getSize());
	standard->getDataAtTime(time, state->getSize(), data);
	for (int j = 0; j < state->getSize(); ++j) {
        stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard->getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j] << endl;
        ASSERT_EQUAL(data[j], state->getData()[j], 1.0e-3, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str();
	}
}

void testGait2354WithController() {
    ForwardTool forward("subject01_Setup_Forward_Controller.xml");
	forward.run();
	Storage results("ResultsCorrectionController/subject01_states.sto");
	Storage* standard = new Storage();
    string statesFileName("std_subject01_walk1_states.sto");
    forward.loadStatesStorage( statesFileName, standard );

	Array<double> data;
	int i = results.getSize() - 1;
	StateVector* state = results.getStateVector(i);
	double time = state->getTime();
	data.setSize(state->getSize());
	standard->getDataAtTime(time, state->getSize(), data);
	int nc = forward.getModel().getNumCoordinates();
	for (int j = 0; j < nc; ++j) {      
	    stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard->getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j];
        ASSERT_EQUAL(data[j], state->getData()[j], 1e-2, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str() << endl;
	}
}

void testPendulumExternalLoadWithPointInGround() {
	ForwardTool forward("pendulum_ext_point_in_ground_Setup_Forward.xml");cout << endl;
	forward.run();
	Storage results("Results/pendulum_ext_gravity_point_in_ground_states.sto");
	ASSERT(results.getFirstTime() == 0.0);
	ASSERT(results.getLastTime() == 1.0);

	Storage standard("Results/pendulum_states.sto");
 

	Array<double> data;
	int i = results.getSize() - 1;
	StateVector* state = results.getStateVector(i);
	double time = state->getTime();
	data.setSize(state->getSize());
	standard.getDataAtTime(time, state->getSize(), data);
	int nc = forward.getModel().getNumCoordinates();
	for (int j = 0; j < nc; ++j) {      
	    stringstream message;
		message << "t=" << time <<" state# "<< j << " " << standard.getColumnLabels()[j+1] << " std=" << data[j] <<"  computed=" << state->getData()[j];
        ASSERT_EQUAL(data[j], state->getData()[j], 1e-2, __FILE__, __LINE__, "ASSERT_EQUAL FAILED " + message.str());
		cout << "ASSERT_EQUAL PASSED " << message.str() << endl;
	}
}

// -------------------------------------------------------------------------------------------------
//  OperatorList::display(StateVector sv)
//
//  Use: Output a text visualisation of the operator list
//  Input: StateVector sv - the spin configuration
//  Output: none
// -------------------------------------------------------------------------------------------------
void OperatorList::display(StateVector sv)
{
    int N = sv.getSize();
    int L = getSize();

    const char* spins[2] = {"U", "D"};
    const char* ops[3] = {" ", "-", "="};

    for (int i = 0; i < L; i++) {
        cout << " ";
        for (int j = 0; j < N; j++)
            cout << spins[sv(j)] << " ";
        cout  << " (" << i << ")" << endl;

        int op = list->operators[i].op;
        if (op == 0) 
            cout << endl;
        else {
            int spin1 = list->operators[i].spin;
            int spin2 = spin1 + 1;
            if (spin2 > (N - 1))
                spin2 = 0;
            
            if (spin2 == 0) {
                cout << " " << ops[op];
                for (int k = 0; k < (2*N - 3); k++)
                    cout << " ";
                cout << ops[op] << endl;
            } else {
                for (int k = 0; k < (2*spin1 + 1); k++)
                    cout << " ";            
                for (int k = 0; k < 3; k++)
                    cout << ops[op];
                cout << endl;
            }

            // flip spins if off diagonal operator
            if (op == 2) {
                sv.flip(spin1);
                sv.flip(spin2);
            }            
        }
    }
    cout << " ";
    for (int j = 0; j < N; j++)
        cout << spins[sv(j)] << " ";
    cout << endl << endl << endl;
}

// -------------------------------------------------------------------------------------------------
//  OperatorList::loopUpdate(StateVector& sv, double beta)
//
//  Use: Performs a loop cluster update on the operator list
//  Input: StateVector& sv - the spin configuration
//         double beta - the inverse temperature
//  Output: none
// -------------------------------------------------------------------------------------------------
void OperatorList::loopUpdate(StateVector& sv, double beta)
{
    int N = sv.getSize();     // Number of spins
    int L = getSize();        // Operator list length

    int noDOD = countDOD();
    if (noDOD > 0) {
        // j is the operator loop index
        int j;
        // Select a random non-unit operator...
        do {
            j = (int)(ran()*L);
            if (j > (L - 1)) 
              j = L - 1;
        } while (list->operators[j].op == 0); 

        // i is the state loop index
        int i = j;

        // direction = -1 : backward direction
        // direction = 1 : forward direction
        int direction = 1;
        int loops = 0;
        int spin = list->operators[j].spin;
        int startop = j;
        int startstate = i;
        int startspin = spin;
        int spin1 = spin;
        int spin2 = (spin1 + 1) % N;

        // change this operator from diagonal to off-diagonal and vice versa...
        changeOp(j);

        // this is the actual loop update
        do {
            if (loops > 0) {
                // Break the loop when we return to the starting operator
                if ((i == startstate) &&(j == startop))
                    if (spin == startspin) {
                        changeOp(j);
                        break;
                    } else if ((spin - 1) == startspin)
                        break;
            }

            j += direction;
            // we use periodic boundary conditions
            if (j < 0) {
                j = L - 1;
                sv.flip(spin);
            }
            if (j > (L - 1)) {
                j = 0;
                sv.flip(spin);
            }

            int op = list->operators[j].op;
            spin1 = list->operators[j].spin;
            spin2 = (spin1 + 1) % N;

            if (op == 0) {
                // in case of a unit operator just update state index
                i += direction;
                if (i < 0) i = L - 1;
                if (i > (L-1)) i = 0;
            } else {
                // diagonal or off-diagonal operator
                if ((spin != spin1) && (spin != spin2)) {
                    i += direction;
                    if (i < 0) i = L - 1;
                    if (i > (L-1)) i = 0;
                } else if (spin == spin1) {
                    spin++;
                    if (spin > (N - 1))
                        spin = 0;
                    direction *= -1;
                    changeOp(j);
                } else if (spin == spin2) {
                    spin--;
                    if (spin < 0)
                        spin = N - 1;
                    direction *= -1;
                    changeOp(j);
                }
            }

            loops = 1;
        } while(true);
    }
}