docs: def LMA(XY,ParIni):
      """
      input: list of x and y values [[x_1, y_1], [x_2, y_2], ....], and a tuple containing an initial guess (a, b, r)
             which is acquired by using an algebraic circle fit (TaubinSVD)
      output: a, b, r.  a and b are the center of the fitting circle, and r is the radius
      %     Geometric circle fit (minimizing orthogonal distances)  
      %     based on the Levenberg-Marquardt scheme in the
      %     "algebraic parameters" A,B,C,D  with constraint B*B+C*C-4*A*D=1
      %        N. Chernov and C. Lesort, "Least squares fitting of circles",
      %        J. Math. Imag. Vision, Vol. 23, 239-251 (2005)
      """
      factorUp=10
      factorDown=0.04
      lambda0=0.01
      epsilon=0.000001
      IterMAX = 50
      AdjustMax = 20
      Xshift=0  
      Yshift=0  
      dX=1  
      dY=0;                                                                                    
      n = len(XY);      # number of data points
  
      anew = ParIni[0] + Xshift
      bnew = ParIni[1] + Yshift
      Anew = old_div(1.,(2.*ParIni[2]))                                                                              
      aabb = anew*anew + bnew*bnew    
      Fnew = (aabb - ParIni[2]*ParIni[2])*Anew 
      Tnew = numpy.arccos(old_div(-anew,numpy.sqrt(aabb))) 
      if bnew > 0:
          Tnew = 2*numpy.pi - Tnew
      VarNew = VarCircle(XY,ParIni) 
  
      VarLambda = lambda0;  
      finish = 0;  
                                                                                                        
      for it in range(0,IterMAX):
                                                                        
          Aold = Anew  
          Fold = Fnew
          Told = Tnew
          VarOld = VarNew
  
          H = numpy.sqrt(1+4*Aold*Fold);                                                                 
          aold = -H*numpy.cos(Told)/(Aold+Aold) - Xshift;
          bold = -H*numpy.sin(Told)/(Aold+Aold) - Yshift;
          Rold = old_div(1,abs(Aold+Aold)); 
  
          DD = 1 + 4*Aold*Fold; 
          D = numpy.sqrt(DD);  
          CT = numpy.cos(Told); 
          ST = numpy.sin(Told);    
          H11=0; 
          H12=0; 
          H13=0; 
          H22=0; 
          H23=0; 
          H33=0; 
          F1=0; 
          F2=0; 
          F3=0;           
                                                              
          for i in range(0,n):
              Xi = XY[i,0] + Xshift;   
              Yi = XY[i,1] + Yshift;       
              Zi = Xi*Xi + Yi*Yi;  
              Ui = Xi*CT + Yi*ST;             
              Vi =-Xi*ST + Yi*CT;
  
              ADF = Aold*Zi + D*Ui + Fold;    
              SQ = numpy.sqrt(4*Aold*ADF + 1);           
              DEN = SQ + 1;                                          
              Gi = 2*ADF/DEN;   
              FACT = 2/DEN*(1 - Aold*Gi/SQ);      
              DGDAi = FACT*(Zi + 2*Fold*Ui/D) - Gi*Gi/SQ;                
              DGDFi = FACT*(2*Aold*Ui/D + 1);
              DGDTi = FACT*D*Vi;    
                                                            
              H11 = H11 + DGDAi*DGDAi;                 
              H12 = H12 + DGDAi*DGDFi;                           
              H13 = H13 + DGDAi*DGDTi;                                          
              H22 = H22 + DGDFi*DGDFi;
              H23 = H23 + DGDFi*DGDTi;                                     
              H33 = H33 + DGDTi*DGDTi;                        
                                                   
              F1 = F1 + Gi*DGDAi; 
              F2 = F2 + Gi*DGDFi;    
              F3 = F3 + Gi*DGDTi;
  
  
          for adjust in range(1,AdjustMax):
                                                
  #             Cholesly decomposition                                     
                                                                         
              G11 = numpy.sqrt(H11 + VarLambda);
              G12 = old_div(H12,G11)                                                                              
              G13 = old_div(H13,G11)
              G22 = numpy.sqrt(H22 + VarLambda - G12*G12);                                                              
              G23 = old_div((H23 - G12*G13),G22);                                             
              G33 = numpy.sqrt(H33 + VarLambda - G13*G13 - G23*G23);                
                                                                                     
              D1 = old_div(F1,G11);                                            
              D2 = old_div((F2 - G12*D1),G22);                                                              
              D3 = old_div((F3 - G13*D1 - G23*D2),G33);                
  
              dT = old_div(D3,G33);  
              dF = old_div((D2 - G23*dT),G22) 
              dA = old_div((D1 - G12*dF - G13*dT),G11) 
                                                                                     
  #            updating the parameters
                                                                                              
              Anew = Aold - dA;  
              Fnew = Fold - dF;                             
              Tnew = Told - dT;
  
              if 1+4*Anew*Fnew < epsilon and VarLambda>1:  
                  Xshift = Xshift + dX;                                          
                  Yshift = Yshift + dY;                                                                               
                                                                                       
                  H = numpy.sqrt(1+4*Aold*Fold);                               
                  aTemp = -H*numpy.cos(Told)/(Aold+Aold) + dX;                                     
                  bTemp = -H*numpy.sin(Told)/(Aold+Aold) + dY;                                      
                  rTemp = old_div(1,abs(Aold+Aold));                                       
                                                                               
                  Anew = old_div(1,(rTemp + rTemp));                         
                  aabb = aTemp*aTemp + bTemp*bTemp;                          
                  Fnew = (aabb - rTemp*rTemp)*Anew;                             
                  Tnew = numpy.arccos(old_div(-aTemp,numpy.sqrt(aabb)));                                       
                  if bTemp > 0:
                      Tnew = 2*numpy.pi - Tnew;           
                  VarNew = VarOld;                                         
                  break;                               
  
              
              if 1+4*Anew*Fnew < epsilon:  
                  VarLambda = VarLambda * factorUp;             
                  continue;              
  
              DD = 1 + 4*Anew*Fnew;                  
              D = numpy.sqrt(DD);                                                         
              CT = numpy.cos(Tnew);                                
              ST = numpy.sin(Tnew);    
                      
              GG = 0;                
                              
  
              for i in range(0, n):
                  Xi = XY[i,0] + Xshift;          
                  Yi = XY[i,1] + Yshift;    
                  Zi = Xi*Xi + Yi*Yi; 
                  Ui = Xi*CT + Yi*ST;            
                                                   
                  ADF = Anew*Zi + D*Ui + Fnew;                    
                  SQ = numpy.sqrt(4*Anew*ADF + 1);               
                  DEN = SQ + 1;                   
                  Gi = 2*ADF/DEN; 
                  GG = GG + Gi*Gi;
                                     
              VarNew = old_div(GG,(n-3));    
           
              H = numpy.sqrt(1+4*Anew*Fnew);               
              anew = -H*numpy.cos(Tnew)/(Anew+Anew) - Xshift;  
              bnew = -H*numpy.sin(Tnew)/(Anew+Anew) - Yshift;  
              Rnew = old_div(1,abs(Anew+Anew)); 
  
              if VarNew <= VarOld: 
                  progress = old_div((abs(anew-aold) + abs(bnew-bold) + abs(Rnew-Rold)),(Rnew+Rold));      
                  if progress < epsilon: 
                      Aold = Anew;          
                      Fold = Fnew;      
                      Told = Tnew;           
                      VarOld = VarNew # %#ok<NASGU>  
                      finish = 1;     
                      break;  
  
                  VarLambda = VarLambda * factorDown
                  break  
              else:                 #    %   no improvement  
                  VarLambda = VarLambda * factorUp;      
                  continue;     
  
          if finish == 1:
              break
  
      H = numpy.sqrt(1+4*Aold*Fold);                                                                                        
      result_a = -H*numpy.cos(Told)/(Aold+Aold) - Xshift;                                                      
      result_b = -H*numpy.sin(Told)/(Aold+Aold) - Yshift;                                                 
      result_r = old_div(1,abs(Aold+Aold));       
  
      return result_a, result_b, result_r