spidx.hpp

来自「Gaussian Mixture Algorithm」· HPP 代码 · 共 626 行 · 第 1/2 页

		} else { if (ndimset) { throw(-8); } }
		if (i6>=0) {
			if (!ndimset) { dim = 7; ndimset = true; }
		} else { if (ndimset) { throw(-7); } }
		if (i5>=0) {
			if (!ndimset) { dim = 6; ndimset = true; }
		} else { if (ndimset) { throw(-6); } }
		if (i4>=0) {
			if (!ndimset) { dim = 5; ndimset = true; }
		} else { if (ndimset) { throw(-5); } }
		if (i3>=0) {
			if (!ndimset) {dim = 4; ndimset = true; }
		} else { if (ndimset) { throw(-4); } }
		if (i2>=0) {
			if (!ndimset) { dim = 3; ndimset = true; }
		} else { if (ndimset) { throw(-3); } }
		if (i1>=0) {
			if (!ndimset) { dim = 2; ndimset = true; }
		} else { if (ndimset) { throw(-2); } }
		if (i0>=0) {
			if (!ndimset) { dim = 1; ndimset = true; }
		} else { if (ndimset) { throw(-1); } }
		if (!ndimset) { throw(0); }
	}
	catch(int v) {return false; }
	return (dim == this->ndim);
}

template<class T> bool spIdx<T>::is_within_bounds(intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7){
	intg bla[8] = { i0, i1, i2, i3, i4, i5, i6, i7};
	for(int i = 0; i < ndim; i++) if((bla[i] < 0) || (bla[i]>=mydims[i])) return false;
	return true;
}

template<class T> intg spIdx<T>::pos_to_index(intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7){
	if(nelem == 0) return -1;
	intg pos = -1;
	bool found = false;
	{ idx_bloop1(ind, *myindex, intg){
		pos++;
		if(pos >= nelem) break;
		intg s0 = ind.get(0);
		if(s0==i0){
			if(ndim >= 2) {
				intg s1 = ind.get(1);
				if(s1==i1){
					if(ndim >= 3) {
						intg s2 = ind.get(2);
						if(s2==i2){
							if(ndim >= 4) {
								intg s3 = ind.get(3);
								if(s3==i3){
									if(ndim >= 5) {
										intg s4 = ind.get(4);
										if(s4==i4){
											if(ndim >= 6) {
												intg s5 = ind.get(5);
												if(s5==i5){
													if(ndim >= 7) {
														intg s6 = ind.get(6);
														if(s6==i6){
															if(ndim >= 8) {
																intg s7 = ind.get(7);
																if(s7==i7){
																	found = true;
																	break;
																}
															} else {
																found = true;
																break;
															}
														}
													} else {
														found = true;
														break;
													}
												}
											} else {
												found = true;
												break;
											}
										}
									} else {
										found = true;
										break;
									}
								}
							} else {
								found = true;
								break;
							}
						}
					} else {
						found = true;
						break;
					}
				}
			} else {
				found = true;
				break;
			}
		}
	}}
	if(!found || (pos == -1)) return -1;
	else return pos;
}

template<class T> void spIdx<T>::index_to_pos(intg pos, intg &i0, intg &i1, intg &i2, intg &i3, intg &i4, intg &i5, intg &i6, intg &i7){
	if((pos<0) || (pos >= nelem)){i0=i1=i2=i3=i4=i5=i6=i7=-1; return;}
	intg* myptr = myindex->idx_ptr();
	const intg* mymod = myindex->mods();
	i0 = (ndim >= 1)? *(myptr + pos * mymod[0]) : -1;
	i1 = (ndim >= 2)? *(myptr + pos * mymod[0] + mymod[1]) : -1;
	i2 = (ndim >= 3)? *(myptr + pos * mymod[0] + 2* mymod[1]) : -1;
	i3 = (ndim >= 4)? *(myptr + pos * mymod[0] + 3* mymod[1]) : -1;
	i4 = (ndim >= 5)? *(myptr + pos * mymod[0] + 4* mymod[1]) : -1;
	i5 = (ndim >= 6)? *(myptr + pos * mymod[0] + 5* mymod[1]) : -1;
	i6 = (ndim >= 7)? *(myptr + pos * mymod[0] + 6* mymod[1]) : -1;
	i7 = (ndim >= 8)? *(myptr + pos * mymod[0] + 7* mymod[1]) : -1;
}

template<class T> T spIdx<T>::get(intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7){
	if(!has_right_dimension(i0, i1, i2, i3, i4, i5, i6, i7)){ylerror("spIdx: bad number of arguments in get"); return (T)BACKGROUND;}
	if(!is_within_bounds(i0, i1, i2, i3, i4, i5, i6, i7)){
		ylerror("spIdx: bad index in get");
		return (T)BACKGROUND;
	}
	intg pos = pos_to_index(i0, i1, i2, i3, i4, i5, i6, i7);
	if(pos == -1) return (T)BACKGROUND;
	else {
		T* myptr = myvalues->idx_ptr();
		intg mymod = myvalues->mod(0);
		return *(myptr + pos * mymod);
	}
}

template<class T> T spIdx<T>::set(T val, intg i0, intg i1, intg i2, intg i3, intg i4, intg i5, intg i6, intg i7){
	if(!has_right_dimension(i0, i1, i2, i3, i4, i5, i6, i7)) {ylerror("spIdx: bad number of arguments in set"); return (T)BACKGROUND; };
	if(!is_within_bounds(i0, i1, i2, i3, i4, i5, i6, i7)){
		ylerror("spIdx: bad index in set");
		return (T)BACKGROUND;
	}
	intg pos = pos_to_index(i0, i1, i2, i3, i4, i5, i6, i7);
	if(pos != -1){
		if(val != BACKGROUND){ myvalues->set(val, pos); return val;}
		else {
			intg N = nelem;
			if(N == 1) { nelem = 0; return val;}
			else {
				intg* myptr = myindex->idx_ptr();
				const intg* mymod = myindex->mods();
				*(myptr + pos * mymod[0]) = *(myptr + (N-1) * mymod[0]);
				if(ndim>=2) *(myptr + pos * mymod[0] + mymod[1]) = *(myptr + (N-1) * mymod[0] + mymod[1]);
				if(ndim>=3) *(myptr + pos * mymod[0] + 2*mymod[1]) = *(myptr + (N-1) * mymod[0] + 2*mymod[1]);
				if(ndim>=4) *(myptr + pos * mymod[0] + 3*mymod[1]) = *(myptr + (N-1) * mymod[0] + 3*mymod[1]);
				if(ndim>=5) *(myptr + pos * mymod[0] + 4*mymod[1]) = *(myptr + (N-1) * mymod[0] + 4*mymod[1]);
				if(ndim>=6) *(myptr + pos * mymod[0] + 5*mymod[1]) = *(myptr + (N-1) * mymod[0] + 5*mymod[1]);
				if(ndim>=7) *(myptr + pos * mymod[0] + 6*mymod[1]) = *(myptr + (N-1) * mymod[0] + 6*mymod[1]);
				if(ndim>=8) *(myptr + pos * mymod[0] + 7*mymod[1]) = *(myptr + (N-1) * mymod[0] + 7*mymod[1]);
				myvalues->set(myvalues->get(N-1), pos);
				nelem--;
				return val;
			}
		}
	} else {
		if(val == BACKGROUND) return val;
		if(nelem == myindex->dim(0)){
			myindex->resize_chunk(100, nelem+1, ndim);
			myvalues->resize_chunk(100, nelem+1);
		}
		intg* myptr = myindex->idx_ptr();
		const intg* mymod = myindex->mods();
		*(myptr + nelem * mymod[0]) = i0;
		if(ndim>=2) *(myptr + nelem * mymod[0] + mymod[1]) = i1;
		if(ndim>=3) *(myptr + nelem * mymod[0] + 2*mymod[1]) = i2;
		if(ndim>=4) *(myptr + nelem * mymod[0] + 3*mymod[1]) = i3;
		if(ndim>=5) *(myptr + nelem * mymod[0] + 4*mymod[1]) = i4;
		if(ndim>=6) *(myptr + nelem * mymod[0] + 5*mymod[1]) = i5;;
		if(ndim>=7) *(myptr + nelem * mymod[0] + 6*mymod[1]) = i6;
		if(ndim>=8) *(myptr + nelem * mymod[0] + 7*mymod[1]) = i7;
		myvalues->set(val, nelem);
		nelem++;
		return val;
	}
}

template<class T> void spIdx<T>::clean(){
	for(intg i = 0; i < nelem; i++){
		if(myvalues->get(i) == BACKGROUND){
			intg i0, i1, i2, i3, i4, i5, i6, i7;
			index_to_pos(i, i0, i1, i2, i3, i4, i5, i6, i7);
			set(BACKGROUND, i0, i1, i2, i3, i4, i5, i6, i7);
			i--;
		}
	}
	int N = nelem > 0 ? nelem : 1;
	myindex->resize(N, myindex->dim(1));
	myvalues->resize(N);
}

template<class T> void spIdx<T>::pretty(FILE *f){
	fprintf(f, "spIdx: at address %ld \n",(intg)this);
	fprintf(f, "  dimensions : ");
	for(int i = 0; i< ndim; i++) fprintf(f, "%li ",(intg)mydims[i]);
	fprintf(f, "\n Number of non-zero elements : %li\n", nelem);
	fprintf(f, "  Index: at address %ld\n",(intg)myindex);
	myindex->pretty(f);
	fprintf(f, "  Values: at address %ld\n",(intg)myvalues);
	myvalues->pretty(f);
}

template<class T> void spIdx<T>::pretty(std::ostream& out){
	out << "spIdx: at address " << (intg)this << "\n";
	out << "  dimensions : ";
	for(int i = 0; i< ndim; i++) out << (intg)mydims[i] << " ";
	out << "\n Number of non-zero elements : " << nelem;
	out << "\n  Index: at address " << (intg)myindex << "\n";
	myindex->pretty(out);
	out << "  Values: at address " << (intg)myvalues << "\n";
	myvalues->pretty(out);
}

template<class T> void spIdx<T>::printElems(){ printElems(std::cout);}

template<class T> void spIdx<T>::printElems( std::ostream& out ){
	if(nelem == 0) { out << "[[] ]"; return; }
	out << "[";
	for(int i = 0; i < nelem; i++){
		out << "[[";
		for(int j = 0; j< ndim; j++){
			out << (int)(myindex->get(i,j));
			if(j != ndim-1 ) out << ";";
		}
		out << "] " << printElems_impl_cast(myvalues->get(i)) << "]";
	}
	out << "]\n";
}

template<class T> void spIdx<T>::printElems( FILE* f ){
	if(nelem == 0) { fprintf(f, "[[] ]"); return; }
	fprintf(f, "[");
	for(int i = 0; i < nelem; i++){
		fprintf(f, "[[");
		for(int j = 0; j< ndim; j++) {
			fprintf(f, "%i",(int)(myindex->get(i,j)));
			if(j != ndim-1) fprintf(f, ";");
		}
		fprintf(f, "] %li]", printElems_impl_cast(myvalues->get(i)));
	}
	fprintf(f, "]\n");
}

template<class T> void spIdx<T>::sort(){
	if(nelem <= 1) return;
	intg* ptr = myindex->idx_ptr();
	T* ptrval = myvalues->idx_ptr();
	int mod = myindex->mod(0);
	int modval = myvalues->mod(0);
	rec_sort(ptr, mod, ptrval, modval, nelem);
}

template<class T> void spIdx<T>::rec_sort(intg *ptrind, int modind, T *ptrval, int modval, intg nind){
	if(nind <= 1) return;
	intg* myptr1 = ptrind + modind;
	T *myptr2 = ptrval + modval;
	intg *supind = new intg[2*nind];
	intg *infind = new intg[2*nind];
	T *supval = new T[nind];
	T *infval = new T[nind];
	intg *infind2 = infind, *supind2 = supind;
	T *infval2 = infval, *supval2 = supval;
	intg nsupind = 0, ninfind = 0;
	for(int i = 0; i < nind - 1; i++){
		if(superior_to(ptrind, myptr1, modind)){
			copy(myptr1, infind2, modind);
			copy(myptr2, infval2, modval);
			infind2 += modind;
			infval2+= modval;
			ninfind++;
		} else {
			copy(myptr1, supind2, modind);
			copy(myptr2, supval2, modval);
			supind2 += modind;
			supval2 += modval;
			nsupind++;
		}
		myptr1 += modind;
		myptr2 += modval;
	}
	assert(ninfind+nsupind+1 == nind);
	copy( ptrind, ptrind + (ninfind*modind), 1*modind);
	copy( ptrval, ptrval + (ninfind*modval), 1*modval);
	copy(infind, ptrind, ninfind*modind);
	copy(infval, ptrval, ninfind*modval);
	copy(supind, ptrind + modind*(ninfind +1), nsupind*modind);
	copy(supval, ptrval + modval*(ninfind +1), nsupind*modval);
	delete infind;
	delete supind;
	delete infval;
	delete supval;
	rec_sort(ptrind, modind, ptrval, modval, ninfind);
	rec_sort(ptrind+modind*(ninfind+1), modind, ptrval+modval*(ninfind+1), modval, nsupind);
}

template<class T> void copy(T* pt1, T* pt2, int ndim){
	for(int i = 0; i < ndim; i++) pt2[i] = pt1[i];
}


}

#endif /* SPIDX_HPP_ */