FITTING THE INTERNAL LENGTHS BY LEAST SQUARES

This code was generated automatically from a Maple program. The entries in the matrix do not follow a very regular pattern, as the lengths to be fit are sums or arbitrary paths. In symbolic form, this is rather trivial to write. The code includes an iteration due to the fact that negative lengths are not possible. Negative lengths are converted to 0 and the rest of the lengths are re-fitted. The value of the least-squares is computed at the end, completely automatically and optimized.

double FiveSubtree( double Wab, double Mab, double Wac, double Mac,
 double Wad, double Mad, double Wae, double Mae,
 double Wbc, double Mbc, double Wbd, double Mbd,
 double Wbe, double Mbe, double Wcd, double Mcd,
 double Wce, double Mce, double Wde, double Mde,
  double *L)
{ double Z1[49], Z2[49], Z3[49], B1[7], B2[7], B3[7], t;
  int i, it, j;
Z1[0] = Z2[0] = Wab+Wac+Wad+Wae;
Z1[1] = Z2[1] = Wab;
Z1[2] = Z2[2] = Wac+Wad+Wae;
Z1[3] = Z2[3] = Wac;
Z1[4] = Z2[4] = Wad+Wae;
Z1[5] = Z2[5] = Wad;
Z1[6] = Z2[6] = Wae;
B1[0] = B2[0] = Mab+Mac+Mad+Mae;
Z1[7] = Z2[7] = Wab;
Z1[8] = Z2[8] = Wab+Wbc+Wbd+Wbe;
Z1[9] = Z2[9] = Wbc+Wbd+Wbe;
Z1[10] = Z2[10] = Wbc;
Z1[11] = Z2[11] = Wbd+Wbe;
Z1[12] = Z2[12] = Wbd;
Z1[13] = Z2[13] = Wbe;
B1[1] = B2[1] = Mab+Mbc+Mbd+Mbe;
Z1[14] = Z2[14] = Wac+Wad+Wae;
Z1[15] = Z2[15] = Wbc+Wbd+Wbe;
Z1[16] = Z2[16] = Wac+Wad+Wbc+Wbd+Wae+Wbe;
Z1[17] = Z2[17] = Wac+Wbc;
Z1[18] = Z2[18] = Wad+Wbd+Wae+Wbe;
Z1[19] = Z2[19] = Wad+Wbd;
Z1[20] = Z2[20] = Wae+Wbe;
B1[2] = B2[2] = Mac+Mad+Mbc+Mbd+Mae+Mbe;
Z1[21] = Z2[21] = Wac;
Z1[22] = Z2[22] = Wbc;
Z1[23] = Z2[23] = Wac+Wbc;
Z1[24] = Z2[24] = Wac+Wbc+Wcd+Wce;
Z1[25] = Z2[25] = Wcd+Wce;
Z1[26] = Z2[26] = Wcd;
Z1[27] = Z2[27] = Wce;
B1[3] = B2[3] = Mac+Mbc+Mcd+Mce;
Z1[28] = Z2[28] = Wad+Wae;
Z1[29] = Z2[29] = Wbd+Wbe;
Z1[30] = Z2[30] = Wad+Wbd+Wae+Wbe;
Z1[31] = Z2[31] = Wcd+Wce;
Z1[32] = Z2[32] = Wad+Wbd+Wcd+Wae+Wbe+Wce;
Z1[33] = Z2[33] = Wad+Wbd+Wcd;
Z1[34] = Z2[34] = Wae+Wbe+Wce;
B1[4] = B2[4] = Mad+Mbd+Mcd+Mae+Mbe+Mce;
Z1[35] = Z2[35] = Wad;
Z1[36] = Z2[36] = Wbd;
Z1[37] = Z2[37] = Wad+Wbd;
Z1[38] = Z2[38] = Wcd;
Z1[39] = Z2[39] = Wad+Wbd+Wcd;
Z1[40] = Z2[40] = Wad+Wbd+Wcd+Wde;
Z1[41] = Z2[41] = Wde;
B1[5] = B2[5] = Mad+Mbd+Mcd+Mde;
Z1[42] = Z2[42] = Wae;
Z1[43] = Z2[43] = Wbe;
Z1[44] = Z2[44] = Wae+Wbe;
Z1[45] = Z2[45] = Wce;
Z1[46] = Z2[46] = Wae+Wbe+Wce;
Z1[47] = Z2[47] = Wde;
Z1[48] = Z2[48] = Wae+Wbe+Wce+Wde;
B1[6] = B2[6] = Mae+Mbe+Mce+Mde;
for( it=0; it<14; it++ ) {
   for( i=0; i<49; i++ ) Z3[i] = Z2[i];
   for( i=0; i<7; i++ ) B3[i] = B2[i];
   gausselimi(Z3,B3,L,7);
   for( i=0; i<7; i++ ) {
      if( L[i] < MinLen-1e-10 ) {
           for( j=0; j<7; j++ ) Z2[7*i+j] = 0;
           Z2[8*i] = 1; B2[i] = MinLen;
           break; }
      else if( fabs(L[i]-MinLen)<1e-10 ) {
               t = -B1[i];
               for( j=0; j<7; j++ ) t += Z1[7*i+j]*L[j];
               if( t<-1e-10 ) {
                  for( j=0; j<7; j++ ) Z2[7*i+j] = Z1[7*i+j];
                  B2[i] = B1[i];
                  break; }
           }
      }
   if( i>=7 ) break;
   }
{ double t1, t4, t7, t10, t13, t16, t19, t22, t29, t36;

      t1 = L[2]+L[4]+L[0]+L[5];
      t4 = L[4]+L[2]+L[1]+L[5];
      t7 = L[2]+L[4]+L[6]+L[0];
      t10 = L[2]+L[4]+L[6]+L[1];
      t13 = L[2]+L[3]+L[0];
      t16 = L[3]+L[2]+L[1];
      t19 = L[4]+L[3]+L[5];
      t22 = L[6]+L[4]+L[3];
      t29 = L[0]+L[1];
      t36 = L[5]+L[6];
      return( t1*t1*Wad+t4*t4*Wbd+t7*t7*Wae+t10*t10*Wbe+t13*t13*Wac+
	t16*t16*Wbc+t19*t19*Wcd+t22*t22*Wce+t29*t29*Wab+t36*t36*Wde -
	2*(t1*Mad+t4*Mbd+t7*Mae+t10*Mbe+t13*Mac+t16*Mbc+t19*Mcd+t22*Mce+
	t29*Mab+Mde*t36) );
  }
}