Golang Matrix.Dims方法代码示例

//MomentTensor returns the moment tensor for a matrix A of coordinates and a column
//vector mass with the respective massess.
func MomentTensor(A *v3.Matrix, massslice []float64) (*v3.Matrix, error) {
	ar, ac := A.Dims()
	var err error
	var mass *mat64.Dense
	if massslice == nil {
		mass = gnOnes(ar, 1)
	} else {
		mass = mat64.NewDense(ar, 1, massslice)
		//		if err != nil {
		//			return nil, err
		//		}
	}
	center, _, err := MassCenter(A, v3.Dense2Matrix(gnCopy(A)), mass)
	if err != nil {
		return nil, errDecorate(err, "MomentTensor")
	}
	sqrmass := gnZeros(ar, 1)
	//	sqrmass.Pow(mass,0.5)
	pow(mass, sqrmass, 0.5) //the result is stored in sqrmass
	//	fmt.Println(center,sqrmass) ////////////////////////
	center.ScaleByCol(center, sqrmass)
	//	fmt.Println(center,"scaled center")
	centerT := gnZeros(ac, ar)
	centerT.Copy(center.T())
	moment := gnMul(centerT, center)
	return v3.Dense2Matrix(moment), nil
}

//Next Reads the next frame in a XTCObj that has been initialized for read
//With initread. If keep is true, returns a pointer to matrix.DenseMatrix
//With the coordinates read, otherwiser, it discards the coordinates and
//returns nil.
func (X *XTCObj) Next(output *v3.Matrix) error {
	if !X.Readable() {
		return Error{TrajUnIni, X.filename, []string{"Next"}, true}
	}
	cnatoms := C.int(X.natoms)
	worked := C.get_coords(X.fp, &X.cCoords[0], cnatoms)
	if worked == 11 {
		X.readable = false
		return newlastFrameError(X.filename, "Next") //This is not really an error and should be catched in the calling function
	}
	if worked != 0 {
		X.readable = false
		return Error{ReadError, X.filename, []string{"Next"}, true}
	}
	if output != nil { //col the frame
		r, c := output.Dims()
		if r < (X.natoms) {
			panic("Buffer v3.Matrix too small to hold trajectory frame")
		}
		for j := 0; j < r; j++ {
			for k := 0; k < c; k++ {
				l := k + (3 * j)
				output.Set(j, k, (10 * float64(X.cCoords[l]))) //nm to Angstroms
			}
		}
		return nil
	}
	return nil //Just drop the frame
}

//MassCenter centers in in the center of mass of oref. Mass must be
//A column vector. Returns the centered matrix and the displacement matrix.
func MassCenter(in, oref *v3.Matrix, mass *mat64.Dense) (*v3.Matrix, *v3.Matrix, error) {
	or, _ := oref.Dims()
	ir, _ := in.Dims()
	if mass == nil { //just obtain the geometric center
		tmp := ones(or)
		mass = mat64.NewDense(or, 1, tmp) //gnOnes(or, 1)
	}
	ref := v3.Zeros(or)
	ref.Copy(oref)
	gnOnesvector := gnOnes(1, or)
	f := func() { ref.ScaleByCol(ref, mass) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return nil, nil, CError{err.Error(), []string{"v3.Matrix.ScaleByCol", "MassCenter"}}
	}
	ref2 := v3.Zeros(1)
	g := func() { ref2.Mul(gnOnesvector, ref) }
	if err := gnMaybe(gnPanicker(g)); err != nil {
		return nil, nil, CError{err.Error(), []string{"v3.gOnesVector", "MassCenter"}}
	}
	ref2.Scale(1.0/mass.Sum(), ref2)
	returned := v3.Zeros(ir)
	returned.Copy(in)
	returned.SubVec(returned, ref2)
	/*	for i := 0; i < ir; i++ {
			if err := returned.GetRowVector(i).Subtract(ref2); err != nil {
				return nil, nil, err
			}
		}
	*/
	return returned, ref2, nil
}

//BestPlane returns a row vector that is normal to the plane that best contains the molecule
//if passed a nil Masser, it will simply set all masses to 1.
func BestPlane(coords *v3.Matrix, mol Masser) (*v3.Matrix, error) {
	var err error
	var Mmass []float64
	cr, _ := coords.Dims()
	if mol != nil {
		Mmass, err = mol.Masses()
		if err != nil {
			return nil, errDecorate(err, "BestPlane")
		}
		if len(Mmass) != cr {
			return nil, CError{fmt.Sprintf("Inconsistent coordinates(%d)/atoms(%d)", len(Mmass), cr), []string{"BestPlane"}}
		}
	}
	moment, err := MomentTensor(coords, Mmass)
	if err != nil {
		return nil, errDecorate(err, "BestPlane")
	}
	evecs, _, err := v3.EigenWrap(moment, appzero)
	if err != nil {
		return nil, errDecorate(err, "BestPlane")
	}
	normal, err := BestPlaneP(evecs)
	if err != nil {
		return nil, errDecorate(err, "BestPlane")
	}
	//MomentTensor(, mass)
	return normal, err
}

//EulerRotateAbout uses Euler angles to rotate the coordinates in coordsorig around by angle
//radians around the axis given by the vector axis. It returns the rotated coordsorig,
//since the original is not affected. It seems more clunky than the RotateAbout, which uses Clifford algebra.
//I leave it for benchmark, mostly, and might remove it later. There is no test for this function!
func EulerRotateAbout(coordsorig, ax1, ax2 *v3.Matrix, angle float64) (*v3.Matrix, error) {
	r, _ := coordsorig.Dims()
	coords := v3.Zeros(r)
	translation := v3.Zeros(ax1.NVecs())
	translation.Copy(ax1)
	axis := v3.Zeros(ax2.NVecs())
	axis.Sub(ax2, ax1) //now it became the rotation axis
	f := func() { coords.SubVec(coordsorig, translation) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return nil, CError{err.Error(), []string{"v3.Matrix.Subvec", "EulerRotateAbout"}}

	}
	Zswitch := RotatorToNewZ(axis)
	coords.Mul(coords, Zswitch) //rotated
	Zrot, err := RotatorAroundZ(angle)
	if err != nil {
		return nil, errDecorate(err, "EulerRotateAbout")
	}
	//	Zsr, _ := Zswitch.Dims()
	//	RevZ := v3.Zeros(Zsr)
	RevZ, err := gnInverse(Zswitch)
	if err != nil {
		return nil, errDecorate(err, "EulerRotateAbout")
	}
	coords.Mul(coords, Zrot) //rotated
	coords.Mul(coords, RevZ)
	coords.AddVec(coords, translation)
	return coords, nil
}

//PDBStringWrite writes a string in PDB format for a given reference, coordinate set and bfactor set, which must match each other
//returns the written string and error or nil.
func PDBStringWrite(coords *v3.Matrix, mol Atomer, bfact []float64) (string, error) {
	if bfact == nil {
		bfact = make([]float64, mol.Len())
	}
	cr, _ := coords.Dims()
	br := len(bfact)
	if cr != mol.Len() || cr != br {
		return "", CError{"Ref and Coords and/or Bfactors dont have the same number of atoms", []string{"PDBStringWrite"}}
	}
	chainprev := mol.Atom(0).Chain //this is to know when the chain changes.
	var outline string
	var outstring string
	var err error
	for i := 0; i < mol.Len(); i++ {
		//	r,c:=coords.Dims()
		//	fmt.Println("IIIIIIIIIIIi", i,coords,r,c, "lllllll")
		writecoord := coords.VecView(i)
		outline, chainprev, err = writePDBLine(mol.Atom(i), writecoord, bfact[i], chainprev)
		if err != nil {
			return "", errDecorate(err, "PDBStringWrite "+fmt.Sprintf("Could not print PDB line: %d", i))
		}
		outstring = strings.Join([]string{outstring, outline}, "")
	}
	outstring = strings.Join([]string{outstring, "END\n"}, "")
	return outstring, nil
}

func pdbWrite(out io.Writer, coords *v3.Matrix, mol Atomer, bfact []float64) error {
	if bfact == nil {
		bfact = make([]float64, mol.Len())
	}
	cr, _ := coords.Dims()
	br := len(bfact)
	if cr != mol.Len() || cr != br {
		return CError{"Ref and Coords and/or Bfactors dont have the same number of atoms", []string{"pdbWrite"}}
	}
	chainprev := mol.Atom(0).Chain //this is to know when the chain changes.
	var outline string
	var err error
	iowriteError := func(err error) error {
		return CError{"Failed to write in io.Writer" + err.Error(), []string{"io.Write.Write", "pdbWrite"}}
	}
	for i := 0; i < mol.Len(); i++ {
		//	r,c:=coords.Dims()
		//	fmt.Println("IIIIIIIIIIIi", i,coords,r,c, "lllllll")
		writecoord := coords.VecView(i)
		outline, chainprev, err = writePDBLine(mol.Atom(i), writecoord, bfact[i], chainprev)
		if err != nil {
			return errDecorate(err, "pdbWrite "+fmt.Sprintf("Could not print PDB line: %d", i))
		}
		_, err := out.Write([]byte(outline))
		if err != nil {
			return iowriteError(err)
		}
	}
	_, err = out.Write([]byte("END")) //no newline, this is in case the write is part of a PDB and one needs to write "ENDMODEL".
	if err != nil {
		return iowriteError(err)
	}
	return nil
}

//Projection returns the projection of test in ref.
func Projection(test, ref *v3.Matrix) *v3.Matrix {
	rr, _ := ref.Dims()
	Uref := v3.Zeros(rr)
	Uref.Unit(ref)
	scalar := test.Dot(Uref) //math.Abs(la)*math.Cos(angle)
	Uref.Scale(scalar, Uref)
	return Uref
}

//AntiProjection returns a vector in the direction of ref with the magnitude of
//a vector A would have if |test| was the magnitude of its projection
//in the direction of test.
func AntiProjection(test, ref *v3.Matrix) *v3.Matrix {
	rr, _ := ref.Dims()
	testnorm := test.Norm(0)
	Uref := v3.Zeros(rr)
	Uref.Unit(ref)
	scalar := test.Dot(Uref)
	scalar = (testnorm * testnorm) / scalar
	Uref.Scale(scalar, Uref)
	return Uref
}

//BestPlaneP takes sorted evecs, according to the eval,s and returns a row vector that is normal to the
//Plane that best contains the molecule. Notice that the function can't possibly check
//that the vectors are sorted. The P at the end of the name is for Performance. If
//That is not an issue it is safer to use the BestPlane function that wraps this one.
func BestPlaneP(evecs *v3.Matrix) (*v3.Matrix, error) {
	evr, evc := evecs.Dims()
	if evr != 3 || evc != 3 {
		return evecs, CError{"goChem: Eigenvectors matrix must be 3x3", []string{"BestPlaneP"}} //maybe this should be a panic
	}
	v1 := evecs.VecView(2)
	v2 := evecs.VecView(1)
	normal := v3.Zeros(1)
	normal.Cross(v1, v2)
	return normal, nil
}

//RotatorToNewY takes a set of coordinates (mol) and a vector (y). It returns
//a rotation matrix that, when applied to mol, will rotate it around the Z axis
//in such a way that the projection of newy in the XY plane will be aligned with
//the Y axis.
func RotatorAroundZToNewY(newy *v3.Matrix) (*v3.Matrix, error) {
	nr, nc := newy.Dims()
	if nc != 3 || nr != 1 {
		return nil, CError{"Wrong newy vector", []string{"RotatorAroundZtoNewY"}}
	}
	if nc != 3 {
		return nil, CError{"Wrong mol vector", []string{"RotatorAroundZtoNewY"}} //this one doesn't seem reachable

	}
	gamma := math.Atan2(newy.At(0, 0), newy.At(0, 1))
	singamma := math.Sin(gamma)
	cosgamma := math.Cos(gamma)
	operator := []float64{cosgamma, singamma, 0,
		-singamma, cosgamma, 0,
		0, 0, 1}
	return v3.NewMatrix(operator)

}

//CenterOfMass returns the center of mass the atoms represented by the coordinates in geometry
//and the masses in mass, and an error. If mass is nil, it calculates the geometric center
func CenterOfMass(geometry *v3.Matrix, mass *mat64.Dense) (*v3.Matrix, error) {
	if geometry == nil {
		return nil, CError{"goChem: nil matrix to get the center of mass", []string{"CenterOfMass"}}
	}
	gr, _ := geometry.Dims()
	if mass == nil { //just obtain the geometric center
		tmp := ones(gr)
		mass = mat64.NewDense(gr, 1, tmp) //gnOnes(gr, 1)
	}
	tmp2 := ones(gr)
	gnOnesvector := mat64.NewDense(1, gr, tmp2) //gnOnes(1, gr)

	ref := v3.Zeros(gr)
	ref.ScaleByCol(geometry, mass)
	ref2 := v3.Zeros(1)
	ref2.Mul(gnOnesvector, ref)
	ref2.Scale(1.0/mass.Sum(), ref2)
	return ref2, nil
}

//RotatorToNewZ takes a matrix a row vector (newz).
//It returns a linear operator such that, when applied to a matrix mol ( with the operator on the right side)
//it will rotate mol such that the z axis is aligned with newz.
func RotatorToNewZ(newz *v3.Matrix) *v3.Matrix {
	r, c := newz.Dims()
	if c != 3 || r != 1 {
		panic("Wrong newz vector")
	}
	normxy := math.Sqrt(math.Pow(newz.At(0, 0), 2) + math.Pow(newz.At(0, 1), 2))
	theta := math.Atan2(normxy, newz.At(0, 2))      //Around the new y
	phi := math.Atan2(newz.At(0, 1), newz.At(0, 0)) //First around z
	psi := 0.000000000000                           // second around z
	sinphi := math.Sin(phi)
	cosphi := math.Cos(phi)
	sintheta := math.Sin(theta)
	costheta := math.Cos(theta)
	sinpsi := math.Sin(psi)
	cospsi := math.Cos(psi)
	operator := []float64{cosphi*costheta*cospsi - sinphi*sinpsi, -sinphi*cospsi - cosphi*costheta*sinpsi, cosphi * sintheta,
		sinphi*costheta*cospsi + cosphi*sinpsi, -sinphi*costheta*sinpsi + cosphi*cospsi, sintheta * sinphi,
		-sintheta * cospsi, sintheta * sinpsi, costheta}
	finalop, _ := v3.NewMatrix(operator) //we are hardcoding opperator so it must have the right dimensions.
	return finalop

}

// RamaCalc Obtains the values for the phi and psi dihedrals indicated in []Ramaset, for the
// structure M.  It returns a slice of 2-element slices, one for the phi the next for the psi
// dihedral, a and an error or nil.
func RamaCalc(M *v3.Matrix, dihedrals []RamaSet) ([][]float64, error) {
	if M == nil || dihedrals == nil {
		return nil, Error{ErrNilData, "", "RamaCalc", "", true}
	}
	r, _ := M.Dims()
	Rama := make([][]float64, 0, len(dihedrals))
	for _, j := range dihedrals {
		if j.Npost >= r {
			return nil, Error{ErrOutOfRange, "", "RamaCalc", "", true}
		}
		Cprev := M.VecView(j.Cprev)
		N := M.VecView(j.N)
		Ca := M.VecView(j.Ca)
		C := M.VecView(j.C)
		Npost := M.VecView(j.Npost)
		phi := chem.Dihedral(Cprev, N, Ca, C)
		psi := chem.Dihedral(N, Ca, C, Npost)
		temp := []float64{phi * (180 / math.Pi), psi * (180 / math.Pi)}
		Rama = append(Rama, temp)
	}
	return Rama, nil
}

//RMSD returns the RSMD (root of the mean square deviation) for the sets of cartesian
//coordinates in test and template.
func RMSD(test, template *v3.Matrix) (float64, error) {
	//This is a VERY naive implementation.
	tmr, tmc := template.Dims()
	tsr, tsc := test.Dims()
	if tmr != tsr || tmc != 3 || tsc != 3 {
		return 0, fmt.Errorf("Ill formed matrices for RMSD calculation")
	}
	tr := tmr
	ctempla := v3.Zeros(template.NVecs())
	ctempla.Copy(template)
	//the maybe thing might not be needed since we check the dimensions before.
	f := func() { ctempla.Sub(ctempla, test) }
	if err := gnMaybe(gnPanicker(f)); err != nil {
		return 0, CError{err.Error(), []string{"v3.Matrix.Sub", "RMSD"}}
	}
	var RMSD float64
	for i := 0; i < template.NVecs(); i++ {
		temp := ctempla.VecView(i)
		RMSD += math.Pow(temp.Norm(0), 2)
	}
	RMSD = RMSD / float64(tr)
	RMSD = math.Sqrt(RMSD)
	return RMSD, nil
}