C++ HCDT::GetNy方法代码示例

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_evaluate(HCDT dat, HCDT idat, HCDT jdat, HCDT kdat, int norm)
{
	if(!idat || (jdat && jdat->GetNN()!=idat->GetNN()) || (kdat && kdat->GetNN()!=idat->GetNN()))	return 0;
	const mglData *dd=dynamic_cast<const mglData *>(dat);
	const mglDataC *dc=dynamic_cast<const mglDataC *>(dat);
	long nx=dat->GetNx(), ny=dat->GetNy(), nz=dat->GetNz();
	mglDataC *r=new mglDataC(idat->GetNx(),idat->GetNy(),idat->GetNz());
	if(dd)
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<mreal>(dd->a,nx,ny,nz, x,y,z):NAN;
		}
	else if(dc)
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<dual>(dc->a,nx,ny,nz, x,y,z):NAN;
		}
	else
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mgl_data_linear(dat, x,y,z):NAN;;
		}
	return r;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_combine(HCDT d1, HCDT d2)
{
	long n1=d1->GetNy(),n2=d2->GetNx(),nx=d1->GetNx();
	if(d1->GetNz()>1 || (n1>1 && d2->GetNy()>1) || d2->GetNz()>1)	return 0;	// wrong dimensions
	mglDataC *r=new mglDataC;
	bool dim2=true;
	if(n1==1)	{	n1=n2;	n2=d2->GetNy();	dim2 = false;	}
	r->Create(nx,n1,n2);
	if(dim2)	n1*=nx;	else	{	n2*=n1;	n1=nx;	}

	const mglDataC *c1=dynamic_cast<const mglDataC *>(d1);
	const mglDataC *c2=dynamic_cast<const mglDataC *>(d2);
	if(c1 && c2)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = c1->a[i]*c2->a[j];
	else if(c1)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = c1->a[i]*d2->vthr(j);
	else if(c2)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = d1->vthr(i)*c2->a[j];
	else
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = d1->vthr(i)*d2->vthr(j);
	return r;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_column(HCDT dat, const char *eq)
{
	const mglData *dd=dynamic_cast<const mglData *>(dat);
	std::vector<mglDataA*> list;
	if(dd && dd->id.length()>0)	for(size_t i=0;i<dd->id.length();i++)
	{
		mglDataT *col = new mglDataT(*dat);
		col->SetInd(i,dd->id[i]);
		list.push_back(col);
	}
	const mglDataC *dc=dynamic_cast<const mglDataC *>(dat);
	if(dc && dc->id.length()>0)	for(size_t i=0;i<dc->id.length();i++)
	{
		mglDataT *col = new mglDataT(*dat);
		col->SetInd(i,dc->id[i]);
		list.push_back(col);
	}
	if(list.size()==0)	return 0;	// no named columns
	mglDataV *t = new mglDataV(dat->GetNy(),dat->GetNz());
	t->s=L"#$mgl";	list.push_back(t);
	mglDataC *r = new mglDataC;
	r->Set(mglFormulaCalcC(eq,list));
	for(size_t i=0;i<list.size();i++)	delete list[i];
	return r;
}

HADT MGL_EXPORT mgl_datac_sum(HCDT dat, const char *dir)
{
	if(!dir || *dir==0)	return 0;
	long nx=dat->GetNx(),ny=dat->GetNy(),nz=dat->GetNz();
	long p[3]={nx,ny,nz};
	dual *b = new dual[nx*ny*nz];
	dual *c = new dual[nx*ny*nz];

	const mglDataC *d=dynamic_cast<const mglDataC *>(dat);
	if(d)	memcpy(c,d->a,nx*ny*nz*sizeof(dual));
	else
#pragma omp parallel for
		for(long i=0;i<nx*ny*nz;i++)	c[i]=dat->vthr(i);

	if(strchr(dir,'z') && nz>1)
	{
		mglStartThreadC(mgl_sumc_z,0,nx*ny,b,c,0,p);
		memcpy(c,b,nx*ny*sizeof(mreal));	p[2] = 1;
	}
	if(strchr(dir,'y') && ny>1)
	{
		mglStartThreadC(mgl_sumc_y,0,nx*p[2],b,c,0,p);
		memcpy(c,b,nx*p[2]*sizeof(mreal));	p[1] = p[2];	p[2] = 1;
	}
	if(strchr(dir,'x') && nx>1)
	{
		mglStartThreadC(mgl_sumc_x,0,p[1]*p[2],b,c,0,p);
		p[0] = p[1];	p[1] = p[2];	p[2] = 1;
	}
	mglDataC *r=new mglDataC(p[0],p[1],p[2]);
	memcpy(r->a,b,p[0]*p[1]*p[2]*sizeof(dual));
	delete []b;	delete []c;	return r;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_trace(HCDT d)
{
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	const mglDataC *dc = dynamic_cast<const mglDataC *>(d);
	mglDataC *r=new mglDataC(nx);
	if(dc)
	{
		if(ny>=nx && nz>=nx)
#pragma omp parallel for
			for(long i=0;i<nx;i++)	r->a[i] = dc->a[i+nx*(i+ny*i)];
		else if(ny>=nx)
#pragma omp parallel for
			for(long i=0;i<nx;i++)	r->a[i] = dc->a[i+nx*i];
		else
#pragma omp parallel for
			for(long i=0;i<nx;i++)	r->a[i] = dc->a[i];
	}
	else if(ny>=nx && nz>=nx)
#pragma omp parallel for
		for(long i=0;i<nx;i++)	r->a[i] = d->v(i,i,i);
	else if(ny>=nx)
#pragma omp parallel for
		for(long i=0;i<nx;i++)	r->a[i] = d->v(i,i);
	else
#pragma omp parallel for
		for(long i=0;i<nx;i++)	r->a[i] = d->v(i);
	return r;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_cont_z_val(HMGL gr, HCDT v, HCDT a, const char *sch, double sv, const char *opt)
{
	long n=a->GetNx(),m=a->GetNy();
	if(n<2 || m<2)	{	gr->SetWarn(mglWarnLow,"ContZ");	return;	}
	gr->SaveState(opt);
	if(mgl_isnan(sv))	sv = gr->GetOrgZ('z');
	if(sv<gr->Min.z || sv>gr->Max.z)	{	gr->SetWarn(mglWarnSlc,"ContZ");	gr->LoadState();	return;	}
	static int cgid=1;	gr->StartGroup("ContZ",cgid++);
	mglDataV xx,yy,zz;	mglData aa;

	int text=0;
	if(mglchr(sch,'t'))	text=1;
	if(mglchr(sch,'T'))	text=2;
	long ss=gr->AddTexture(sch);
	gr->SetPenPal(sch);

	a = fill_slice_z(gr,sv,a,xx,yy,zz,aa);
#pragma omp parallel for
	for(long i=0;i<v->GetNx();i++)
	{
		mreal v0 = v->v(i);
		mgl_cont_gen(gr,v0,a,&xx,&yy,&zz,gr->GetC(ss,v0),text,0);
	}
	gr->EndGroup();
}

//-----------------------------------------------------------------------------
std::vector<mglSegment> MGL_NO_EXPORT mgl_tri_lines(mreal val, HCDT nums, HCDT a, HCDT x, HCDT y, HCDT z)
{
	long n = x->GetNx(), m = nums->GetNy();
	std::vector<mglSegment> lines;
	for(long i=0;i<m;i++)
	{
		register long k1 = long(nums->v(0,i)+0.1);	if(k1<0 || k1>=n)	continue;
		register long k2 = long(nums->v(1,i)+0.1);	if(k2<0 || k2>=n)	continue;
		register long k3 = long(nums->v(2,i)+0.1);	if(k3<0 || k3>=n)	continue;
		register mreal v1 = a->v(k1), v2 = a->v(k2), v3 = a->v(k3);
		register mreal d1 = mgl_d(val,v1,v2), d2 = mgl_d(val,v1,v3), d3 = mgl_d(val,v2,v3);
		mglSegment line;
		if(d1>=0 && d1<=1 && d2>=0 && d2<=1)
		{
			line.p1 = mglPoint(x->v(k1)*(1-d1)+x->v(k2)*d1, y->v(k1)*(1-d1)+y->v(k2)*d1, z->v(k1)*(1-d1)+z->v(k2)*d1);
			line.p2 = mglPoint(x->v(k1)*(1-d2)+x->v(k3)*d2, y->v(k1)*(1-d2)+y->v(k3)*d2, z->v(k1)*(1-d2)+z->v(k3)*d2);
		}
		else if(d1>=0 && d1<=1 && d3>=0 && d3<=1)
		{
			line.p1 = mglPoint(x->v(k1)*(1-d1)+x->v(k2)*d1, y->v(k1)*(1-d1)+y->v(k2)*d1, z->v(k1)*(1-d1)+z->v(k2)*d1);
			line.p2 = mglPoint(x->v(k2)*(1-d3)+x->v(k3)*d3, y->v(k2)*(1-d3)+y->v(k3)*d3, z->v(k2)*(1-d3)+z->v(k3)*d3);
		}
		else if(d3>=0 && d3<=1 && d2>=0 && d2<=1)
		{
			line.p1 = mglPoint(x->v(k1)*(1-d2)+x->v(k3)*d2, y->v(k1)*(1-d2)+y->v(k3)*d2, z->v(k1)*(1-d2)+z->v(k3)*d2);
			line.p2 = mglPoint(x->v(k2)*(1-d3)+x->v(k3)*d3, y->v(k2)*(1-d3)+y->v(k3)*d3, z->v(k2)*(1-d3)+z->v(k3)*d3);
		}
		if(line.p1!=line.p2)	lines.push_back(line);
	}
	return lines;
}

//-----------------------------------------------------------------------------
//
//	TriContV series
//
//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_tricontv_xyzcv(HMGL gr, HCDT v, HCDT nums, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, const char *opt)
{
	mglDataV zz(x->GetNx(),x->GetNy());
	if(!z)	z = &zz;
	if(mgl_check_trig(gr,nums,x,y,z,a,"TriContV"))	return;

	gr->SaveState(opt);
	static int cgid=1;	gr->StartGroup("TriContV",cgid++);
	bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
	long s=gr->AddTexture(sch);
	gr->SetPenPal(sch);

	for(long k=0;k<v->GetNx();k++)
	{
		mreal v0 = v->v(k);		zz.Fill(fixed ? gr->Min.z : v0);
		mreal dv = (gr->Max.c-gr->Min.c)/8, c = gr->GetC(s,v0);
		if(k>0)	dv = v->v(k-1)-v->v(k);
		else if(k<v->GetNx()-1)	dv = v->v(k)-v->v(k+1);
		if(fixed)	dv=-dv;

		const std::vector<mglSegment> curvs = mgl_get_curvs(gr,mgl_tri_lines(v0,nums,a,x,y,fixed?&zz:z));
		for(size_t i=0;i<curvs.size();i++)
		{
			const std::list<mglPoint> &pp=curvs[i].pp;
			long f2=-1,g2=-1;
			for(std::list<mglPoint>::const_iterator it=pp.begin(); it != pp.end(); ++it)
			{
				mglPoint p=*it,q(p.y,-p.x);
				long f1 = f2;	f2 = gr->AddPnt(p,c,q);	p.z+=dv;
				long g1 = g2;	g2 = gr->AddPnt(p,c,q);
				gr->quad_plot(f1,g1,f2,g2);
			}
		}
	}
}

//-----------------------------------------------------------------------------
//
//	Dots series
//
//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_dots_ca(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT c, HCDT a, const char *sch, const char *opt)
{
	long n = x->GetNN(), d, k=1;
	if(x->GetNz()>1) 	k=3;		else if(x->GetNy()>1)	k=2;

	if(y->GetNN()!=n || z->GetNN()!=n || c->GetNN()!=n || (a && a->GetNN()!=n))
	{	gr->SetWarn(mglWarnDim,"Dots");	return;	}
	gr->SaveState(opt);

	d = gr->MeshNum>0 ? mgl_ipow(gr->MeshNum+1,k) : n;
	d = n>d ? n/d:1;

	static int cgid=1;	gr->StartGroup("Dots",cgid++);
	char mk=gr->SetPenPal(sch);
	long ss=gr->AddTexture(sch);
	if(mk==0)	mk='.';
	gr->Reserve(n);

	for(long i=0;i<n;i+=d)
	{
		mglPoint p = mglPoint(x->vthr(i),y->vthr(i),z->vthr(i));
		long pp = gr->AddPnt(p,gr->GetC(ss,c->vthr(i)),mglPoint(NAN),a?gr->GetA(a->vthr(i)):-1);
		gr->mark_plot(pp, mk);
	}
	gr->EndGroup();
}

HMDT MGL_EXPORT mgl_data_flame_2d(HCDT A, HCDT F, long n, long skip)
{
	if(!A || A->GetNx()<7 || n<1)	return 0;	// incompatible dimensions
	if(!F || F->GetNx()<2 || F->GetNz()!=A->GetNy())	return 0;	// incompatible dimensions
	mreal amax=0;
	for(long i=0; i<A->GetNy(); i++)	amax += A->v(6,i);
	if(amax<=0)	return 0;
	
	mglData *f = new mglData(3,n);
	mreal x = 0, y = 0;
	for(long i=0; i<skip; i++)	mgl_flame_2d_point(A, F, x, y,amax);
	for(long i=0; i<n; i++)
	{
		f->a[3*i+2] = mgl_flame_2d_point(A, F, x, y, amax);	// TODO color information ?!!
		f->a[3*i] = x;	f->a[3*i+1] = y;
	}
	return f;
}

//-----------------------------------------------------------------------------
HMDT MGL_EXPORT mgl_datac_arg(HCDT d)
{
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	mglData *r=new mglData(nx,ny,nz);
	const mglDataC *dd = dynamic_cast<const mglDataC*>(d);
	if(dd)
#pragma omp parallel for
		for(long i=0;i<nx*ny*nz;i++)	r->a[i] = arg(dd->a[i]);
	return r;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_dens_z(HMGL gr, HCDT a, const char *sch, double sv, const char *opt)
{
	long n=a->GetNx(),m=a->GetNy();
	if(n<2 || m<2)	{	gr->SetWarn(mglWarnLow,"DensZ");	return;	}
	gr->SaveState(opt);
	if(mgl_isnan(sv))	sv = gr->GetOrgZ('z');
	if(sv<gr->Min.z || sv>gr->Max.z)	{	gr->SetWarn(mglWarnSlc,"DensZ");	gr->LoadState();	return;	}
	mglDataV xx,yy,zz;	mglData aa;
	a = fill_slice_z(gr,sv,a,xx,yy,zz,aa);
	mgl_surf_gen(gr, &xx,&yy,&zz,a, 0, sch);
}

//-----------------------------------------------------------------------------
//
//	IFS series
//
//-----------------------------------------------------------------------------
void MGL_NO_EXPORT mgl_ifs_2d_point(HCDT A, mreal& x, mreal& y, mreal amax)
{
	long i, n=A->GetNy();
	mreal r = amax*mgl_rnd(), sum_prob = 0, x1;
	for(i=0;i<n;i++)
	{
		sum_prob += A->v(6,i);
		if(r<sum_prob)	break;
	}
	x1= A->v(0,i)*x + A->v(1,i)*y + A->v(4,i);
	y = A->v(2,i)*x + A->v(3,i)*y + A->v(5,i);	x = x1;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_datac_set(HADT d, HCDT a)
{
	if(!a)	return;
	const mglDataC *dd = dynamic_cast<const mglDataC *>(a);	// faster for mglData
	mgl_datac_create(d, a->GetNx(), a->GetNy(), a->GetNz());
	if(dd)	// this one should be much faster
		memcpy(d->a, dd->a, d->nx*d->ny*d->nz*sizeof(dual));
	else	// very inefficient!!!
	{
		register long i,j,k;
		for(k=0;k<d->nz;k++)	for(j=0;j<d->ny;j++)	for(i=0;i<d->nx;i++)
			d->a[i+d->nx*(j+d->ny*k)] = a->v(i,j,k);
	}
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_contf_y_val(HMGL gr, HCDT v, HCDT a, const char *sch, double sv, const char *opt)
{
	long n=a->GetNx(),m=a->GetNy();
	if(n<2 || m<2)	{	gr->SetWarn(mglWarnLow,"ContFY");	return;	}
	gr->SaveState(opt);
	if(mgl_isnan(sv))	sv = gr->GetOrgY('y');
	if(sv<gr->Min.y || sv>gr->Max.y)	{	gr->SetWarn(mglWarnSlc,"ContFY");	gr->LoadState();	return;	}
	static int cgid=1;	gr->StartGroup("ContFY",cgid++);
	mglDataV xx,yy,zz;	mglData aa;
	long ss=gr->AddTexture(sch);

	a = fill_slice_y(gr,sv,a,xx,yy,zz,aa);
#pragma omp parallel for
	for(long i=0;i<v->GetNx()-1;i++)
	{
		mreal v0 = v->v(i);
		mgl_contf_gen(gr,v0,v->v(i+1),a,&xx,&yy,&zz,gr->GetC(ss,v0),0);
	}
	gr->EndGroup();
}