# why this simple code is not giving the right answer under Compile

I have two sets of inputs that are to be fed to a simple function. The inputs are essentially a set of triangles and the normal for each triangle. The function will compute the gradient about the points (shown in blue). The notebook can be downloaded from the following Dropbox Link

https://www.dropbox.com/s/yas32nfccd2dzj4/debug%20code.nb?dl=0

The two sets of triangles are supposed to be similar (only separated by a translation in space). I compute the area gradient about point pt1

 grad1 = Block[{ptTri, normal, cross, target, facept,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, source = pt1},
Do[
ptTri = opentr1[[i]];
normal = normOpentr1[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross), {i, 1, Length@normOpentr1}];

Do[
ptTri = closedtri1[[j]];
normal = normClosedtr1[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross), {j, 1, Length@normClosedtr1}];
0.7*closedS + 1*openS
]

(*{0.110728, 0.0466838, 0.752509}*)


Likewise if I compute the area gradient about the second point pt2 I get the same answer.

grad2 = Block[{ptTri, normal, cross, target, facept,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, source = pt2},
Do[
ptTri = opentr2[[i]];
normal = normOpentr2[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross);, {i, 1, Length@normOpentr2}];
Do[
ptTri = closedtr2[[j]];
normal = normClosedtr2[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross),{j, 1, Length@normClosedtr2}];
0.7*closedS + 1*openS
]

(*{0.110728, 0.0466838, 0.752509}*)


The two answers are quite similar as the difference between grad1 and grad2 is {1.83187*10^-15, -5.55112*10^-16,4.44089*10^-16}

HOWEVER

the same code when under Compile does not give the same answers for the two datasets.

With[{epcc = 0.7, epco = 1.},
Compile[{{point, _Real, 1}, {opentr, _Real, 3}, {normalO, _Real,2}, {closedtr, _Real, 3}, {normalC,
_Real, 2}},
Block[{ptTri, source = point, normal, target, facept, cross,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}, nO = normalO,
nC = normalC, OT = opentr, CT = closedtr},
Do[
ptTri = OT[[i]];
normal = nO[[i]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross);
, {i, 1, Length@nO}];
Do[
ptTri = CT[[j]];
normal = nC[[j]];
cross = If[ptTri[[1]] == source,
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross);
, {j, 1, Length@nC}];

epcc*closedS + epco*openS
], CompilationTarget -> "C"]
];


Though the structure of the code is the same I guess, the outputs are completely different.

 surfaceGradFn2[pt1, opentr1, normOpentr1, closedtri1, normClosedtr1]
(*{-0.0157582, 0.386426, 0.582118} -> this answer is now different but I gave the same inputs as before *)

(* {0.110728, 0.0466838, 0.752509} -> answer is same as before with the same inputs *)


What am I doing wrong here? Why does not the compiled code work while the non-compiled version works? I wish to have the compiled version of the code working.

also, please ignore what the code intends to do. If you may, kindly look at why the compiled output is different.Thanks and grateful for your help.

• This could be because compiled code does not support arbitrary precision. What happens in C if you force a number into a data type that cannot represent that number, is you get a completely different number. – C. E. May 31 '20 at 9:32
• @C.E. thanks. Can you kindly propose a workaround this. Because the builtin version is a bit slow. I have noticed that the compiled code works for some datasets and not others ! I am surprised that if it is arbitrary precision then the compiled version should not work on other datasets as well. Anyways, what do you recommend, how should I change these numbers to something that compile can understand – Ali Hashmi May 31 '20 at 9:46
• Btw I have used arbitrary precision in other codes while using Compile and it seems to be working fine for those pieces of codes – Ali Hashmi May 31 '20 at 9:48
• ok, I take it back. After looking at your numbers it doesn't seem like they are small or large enough to cause a problem. – C. E. May 31 '20 at 9:53
• @C.E. Any workaround. Can I somehow convert these to something which does not give this strange behaviour – Ali Hashmi May 31 '20 at 9:57

the problem was in the comparator operator i.e. the IF statement. I have used the difference and Chop as the fix: If[Chop[ptTri[[1]] - source, 10^-8] == {0., 0., 0.}...] rather than If[ptTri[[1]] == source...]

the code below now works.

ClearAll@surfaceGrad;
With[{epcc = 0.7, epco = 1.},
Compile[{{point, _Real, 1}, {opentr, _Real, 3}, {normalO, _Real,
2}, {closedtr, _Real, 3}, {normalC, _Real, 2}},
Block[{ptTri, source = point, normal, target, facept, cross,
openS = {0., 0., 0.}, closedS = {0., 0., 0.}},
Do[
ptTri = opentr[[i]];
normal = normalO[[i]];
cross = If[Chop[ptTri[[1]] - source, 10^-8] == {0., 0., 0.},
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
openS += (0.5*cross), {i, 1, Length@normalO}];
Do[
ptTri = closedtr[[j]];
normal = normalC[[j]];
cross = If[Chop[ptTri[[1]] - source, 10^-8] == {0., 0., 0.},
{target, facept} = {ptTri[[2]], ptTri[[-1]]};
Cross[normal, facept - target],
{target, facept} = {ptTri[[1]], ptTri[[-1]]};
Cross[normal, target - facept]
];
closedS += (0.5*cross), {j, 1, Length@normalC}];
epcc*closedS + epco*openS
], CompilationTarget -> "C"]
]