Message Boards Message Boards

KKT conditions for ellipse and hyperbolas for this app

Radek Drozd
Posted 11 years ago

how to add conditions for ellipse and hyperbolas to this app?

q = {Off[NMinimize::eit]; Off[NMinimize::lstol];
   DynamicModule[{rng, ptk, szumx, szumy, final}, 
    Panel[Column[{"Theta z przedzia?u od 0 do:", 
       RadioButtonBar[
        Dynamic[rng], {2 Pi -> "2[Pi]", Pi -> "[Pi]", 
         Pi/2 -> "[Pi]/2", Pi/8 -> "[Pi]/8"}], , 
       "Ilo?? generowanych punktów:", 
       RadioButtonBar[
        Dynamic[ptk], {5 -> "5", 20 -> "20", 50 -> "50", 
         100 -> "100"}], , "Warto?ci wektora przesuni?cia", 
       "Warto?? xi:" RadioButtonBar[
         Dynamic[szumx], {-0.1 -> "-0.1", -.5 -> "-0.5", -.9 -> 
           "-0.9", -3 -> "-3"}], 
       "Warto?? yj:" RadioButtonBar[
         Dynamic[szumy], {0.1 -> "0.1", .5 -> "0.5", .9 -> "0.9", 
          3 -> "3"}], 
       Button["Generuj Elipse", 
        While[True,(*losowo wybra? akceptowalne wspó?czynniki*){a, b, 
           c, d, f, g} = RandomReal[{-10, 10}, 6];
         [CapitalDelta] = -c d^2 + 2 b d f - a f^2 - b^2 g + a c g;
         j = -b^2 + a c; i = a + c;
         If[[CapitalDelta] != 0 && j > 0 && [CapitalDelta]/i < 0, 
          Break[]]];
        ellipse = a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g;
        (*?rodek elipsy w ogólnej formie {(c d-b f)/(b^2-a c),(a f-
        b d)/(b^2-a c)}*)points = Table[theta = RandomReal[{0, rng}];

          ksol = FindRoot[(ellipse /. {x -> 
                 k*Cos[theta] + (c d - b f)/(b^2 - a c), 
                y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)}) == 
             0, {k, 1.}];

          Point[{x, y}] /. {x -> 
              k*Cos[theta] + (c d - b f)/(b^2 - a c), 
             y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)} /. 
           ksol, {ptk}];
        nearpoints = 
         points /. 
          Point[{x_, y_}] :> 
           Point[{x + RandomReal[{szumx, szumy}], 
             y + RandomReal[{szumx, szumy}]}];
        (*minimalne i maksymalne warto?ci x i y dla elipsy*)
        yplotrange = 
         Flatten[{y, 
           Sort[{(2*b*d - 2*a*f + 
                Sqrt[(2*b*d - 2*a*f)^2 - 
                  4*(b^2 - a*c)*(d^2 - a*g)])/(2*(-b^2 + a*c)), (-2*b*
                 d + 2*a*f + 
                Sqrt[(2*b*d - 2*a*f)^2 - 
                  4*(b^2 - a*c)*(d^2 - a*g)])/(2*(b^2 - a*c))}]}];
        xplotrange = 
         Flatten[{x, 
           Sort[{(2*c*d - 2*b*f + 
                Sqrt[(-2*c*d + 2*b*f)^2 - 
                  4*(b^2 - a*c)*(f^2 - c*g)])/(2*(b^2 - a*c)), (-2*c*
                 d + 2*b*f + 
                Sqrt[(-2*c*d + 2*b*f)^2 - 
                  4*(b^2 - a*c)*(f^2 - c*g)])/(2*(-b^2 + a*c))}]}];
        (*zminimalizowa? dystans pobliskich punktów dla nowej ogólnej \n    elpisy*)nearCoords = nearpoints[[All, 1]];
        {xs, ys} = Transpose[nearCoords];
        newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
        distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
        {res, coes} = NMinimize[distance, {aa, bb, cc, dd, ff, gg}];
        scaleup = 
         FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
        esolve = Expand[scaleup*(newellipse /. coes)];
        final = 
         Show[ContourPlot[{ellipse == 0, esolve == 0}, 
           Evaluate[xplotrange], Evaluate[yplotrange]], 
          Graphics[points], Graphics[{Red, nearpoints}], 
          ImageSize -> {500, Automatic}];], Dynamic@final
       (*,Dynamic@Grid[Join[{{"points","nearby points"}},
       Transpose[{points,nearpoints}][[All,All,1]]],Frame->All,
       Alignment->Left]*)}]]]};
Grid[{{Labeled["(A)", q], Labeled["(B)", q]}, {Labeled["(C)", q], 
   Labeled["(D)", q]}}]
POSTED BY: Radek Drozd
19 Replies
Sort By:

[Long sigh.]
POSTED BY: Daniel Lichtblau

This is still not a clear question. For one, it has a large block of code that is difficult to follow. For another, it is not at all clear what intervention you are seeking for that code.

I will infer that you wish to constraint some best fit to a quadric to give an ellipse rather than a hyperbola. If so, provide a very simple example of inputs, and the code that is failing to produce what you require. Throwing a jumble of code and having readers sort through to find an example of some vaguely specified misbehavior is just not a reasonable thing to do. Most especially given that it does not cut and paste to something that will parse.

If you cannot provide a basic example, then likely this thread will get closed as being outside the scope of this forum.
POSTED BY: Daniel Lichtblau

Let me see if I follow this. You have a set of points that you wish to fit to an ellipse. And you are trying to avoid accidently fitting to a hyperbola. Is that it? If so, provide the following.

A sample set of such points, where you get a bad result (a hyperbola).

The specific code you use to fit this point set to an ellipse.

The undesired outcome.
POSTED BY: Daniel Lichtblau
Bill Simpson
Bill Simpson
Posted 11 years ago

Perhaps some generous, patient and skilled individual could engage the student in a protracted intense Socratic dialog and thus enable him to discover and understand what 67 replies to and 11934 views of "help me" have not been able to accomplish. I fear the responses thus far are training the student to be bad instead of good.
POSTED BY: Bill Simpson

And an upvote goes to Bill Simpson. For an elliptic sense of hyperbole, I think.
POSTED BY: Daniel Lichtblau

You need a mathematical condition that distinguishes an ellipse from hyperbola. Requiring that the curve only goes a finite distance from the origin or its center is one possibility.
POSTED BY: Frank Kampas

Before my hair catches fire, have a look at

http://www.wiley.com/college/anton/0471482749/appendices/app_h.pdf

In brief: discriminant<0 is the ellipse condition. Adding disc<=epsilon to the NMinimize, for some modest positive epsilon, should suffice to give only ellipses. It's just more than I'm willing to check at this point.
POSTED BY: Daniel Lichtblau

I thought there must be some simple criterion.

In[1]:= data = Table[{x, 1/x}, {x, Range[6]}]

Out[1]= {{1, 1}, {2, 1/2}, {3, 1/3}, {4, 1/4}, {5, 1/5}, {6, 1/6}}

In[5]:= val = (a x^2 + b*x*y + c*y^2 + d*x + e*y + f)^2 /. 
  Table[Thread[{x, y} -> data[[i]]], {i, Length[data]}]

Out[5]= {(a + b + c + d + e + f)^2, (4 a + b + c/4 + 2 d + e/2 + 
   f)^2, (9 a + b + c/9 + 3 d + e/3 + f)^2, (16 a + b + c/16 + 4 d + 
   e/4 + f)^2, (25 a + b + c/25 + 5 d + e/5 + f)^2, (36 a + b + c/
   36 + 6 d + e/6 + f)^2}

In[10]:= nmres = 
 NMinimize[{Plus @@ val, b^2 - 4*a*c <= -.01}, {a, b, c, d, e, f}]

Out[10]= {0.0000264096, {a -> 0.0071953, b -> 2.4893*10^-19, 
  c -> 0.347449, d -> -0.0922589, e -> -0.658444, f -> 0.396278}}

Show[ContourPlot[
  Evaluate[a x^2 + b*x*y + c y^2 + d*x + e*y + f /. nmres[[2]]] == 
   0, {x, 0, 16}, {y, -4, 4}], ListPlot[data]]

gives an ellipse that passes through the data points.
POSTED BY: Frank Kampas
Radek Drozd
Radek Drozd
Posted 11 years ago

this app application: - Generate elipse - Generate points on the ellipse - Add noise to the ellipse - Adapts elipse of noisy points (optimize) My question is: - I mean KKT adding terms to always go out an ellipse or a second condition - the hyperbole
POSTED BY: Radek Drozd
Radek Drozd
Radek Drozd
Posted 11 years ago

yes i want to fit elipse in the points and want to add KKT conditions to always go out ellipse, let's skip the hyperbole
POSTED BY: Radek Drozd
Radek Drozd
Radek Drozd
Posted 11 years ago

somebody will help me with this?
POSTED BY: Radek Drozd

The KKT conditions are for constrained minimization. Isn't it sufficient to minimize the distance from the points to the ellipse? What are the constraints?
POSTED BY: Frank Kampas
Radek Drozd
Radek Drozd
Posted 11 years ago

The only restriction is added quadratic formula and when generate curves from the clouds of points depends on the settings of the application sometimes to not go ellipse, wants to add KKT to get always ellipse,
POSTED BY: Radek Drozd

If I generate 6 points on an hyperbola

data = Table[{x, 1/x}, {x, Range[6]}]

and use NMinimize to minimize the total of 

(a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g)^2

I get

{3.85186*10^-33, {a -> -1.80964*10^-17, b -> -0.116826, 
  c -> -5.97885*10^-16, d -> 1.02054*10^-16, f -> 5.8286*10^-16, 
  g -> 0.233652}}

which corresponds to x*y = 1 as expected. So, is the question, what constraints, on a,b,c,d,f,g force the curve to be an ellipse?
POSTED BY: Frank Kampas
Radek Drozd
Radek Drozd
Posted 11 years ago

why hyperbola? I have the general quadratic formula 

  ellipse = a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g;

with this conditions for the ellipse

  [CapitalDelta] = -c d^2 + 2 b d f - a f^2 - b^2 g + a c g;
         j = -b^2 + a c; i = a + c;
         If[[CapitalDelta] != 0 && j > 0 && [CapitalDelta]/i < 0, Break[]]];

and next i generate points on the generated ellipse

 Point[{x, y}] /. {x -> k*Cos[theta] + (c d - b f)/(b^2 - a c),  y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)} /. ksol, {ptk}];

and i adds noise to the points in order to generate the next points

 nearpoints =  points /.  Point[{x_, y_}] :> Point[{x + RandomReal[{szumx, szumy}], y + RandomReal[{szumx, szumy}]}];

and from nearpoints optimized to create a new elipse but depending on the settings of the program does not always go ellipse because the noise will be too big or too little will be generated points or points spreading to half or a quarter of an ellipse

 nearCoords = nearpoints[[All, 1]];{xs, ys} = Transpose[nearCoords];
        newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
        distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
        {res, coes} = NMinimize[distance, {aa, bb, cc, dd, ff, gg}];
        scaleup = 
         FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
        esolve = Expand[scaleup*(newellipse /. coes)];

Here we minimize the distance to create a curve, and here I have to add KKT to go out always ellipse
POSTED BY: Radek Drozd

So you want the maximum distance from the origin to any point on the fitted curve to be finite?
POSTED BY: Frank Kampas
Radek Drozd
Radek Drozd
Posted 11 years ago

i only want to add KKT conditions to fit alaways ellipse
POSTED BY: Radek Drozd

enter image description here
POSTED BY: Frank Kampas
Radek Drozd
Radek Drozd
Posted 11 years ago

and as if to add it as a condition to create only ellipse 

[CapitalDelta] = -c d^2 + 2 b d f - a f^2 - b^2 g + a c g;
         j = -b^2 + a c; i = a + c;
         If[[CapitalDelta] != 0 && j > 0 && [CapitalDelta]/i < 0, Break[]]];

and if it was of greater amounts of points and actually cloud of points not lying by the ellipse equation

Place a points with noise to the function E(ax^2 + 2bxy + cy^2 + 2dx + 2f*y + g)=0 like (x0,y0)^2, (x1,y1)^2... and summed them

Further reduce the cost of failure with restrictions

?!=0 J>0 ?/I<0

Normalization a^2+b^2+c^2+d^2+f^2+g^2=1 Find an elipse from the cloud of points with KKT
POSTED BY: Radek Drozd

Group Abstract Group Abstract