<commit_before>/*
  teem: <NAME>'s research software
  Copyright (C) 2002, 2001, 2000, 1999, 1998 University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "gage.h"
#include "privateGage.h"

void
_gageSclAnswer (gageContext *ctx, gagePerVolume *pvl) {
  char me[]="_gageSclAnswer";
  unsigned int query;
  gage_t *ans, gmag=0, *hess, *norm, *gvec, *gten, *k1, *k2, sHess[9], curv=0;
  double tmpMat[9], tmpVec[3], hevec[9], heval[3];
  int *offset;

  gage_t len, gp1[3], gp2[3], nPerp[9], nProj[9], ncTen[9];
  double T, N, D;

  query = pvl->query;
  ans = pvl->ans;
  /* convenience pointers for work below */
  offset = gageKindScl->ansOffset;
  hess = ans + offset[gageSclHessian];
  gvec = ans + offset[gageSclGradVec];
  norm = ans + offset[gageSclNormal];
  gten = ans + offset[gageSclGeomTens];
  k1 = ans + offset[gageSclK1];
  k2 = ans + offset[gageSclK2];
  
  if (1 & (query >> gageSclValue)) {
    /* done if doV */
    if (ctx->verbose) {
      fprintf(stderr, "val = % 15.7f\n", (double)(ans[offset[gageSclValue]]));
    }
  }
  if (1 & (query >> gageSclGradVec)) {
    /* done if doD1 */
    if (ctx->verbose) {
      fprintf(stderr, "gvec = ");
      ell3vPRINT(stderr, gvec);
    }
  }
  if (1 & (query >> gageSclGradMag)) {
    /* this is the true value of gradient magnitude */
    gmag = ans[offset[gageSclGradMag]] = sqrt(ELL_3V_DOT(gvec, gvec));
  }

  /* NB: it would seem that gageParmGradMagMin is completely ignored ... */

  if (1 & (query >> gageSclNormal)) {
    if (gmag) {
      ELL_3V_SCALE(norm, 1.0/gmag, gvec);
      /* polishing ... 
      len = sqrt(ELL_3V_DOT(norm, norm));
      ELL_3V_SCALE(norm, 1.0/len, norm);
      */
    } else {
      ELL_3V_COPY(norm, gageZeroNormal);
    }
  }
  if (1 & (query >> gageSclHessian)) {
    /* done if doD2 */
    if (ctx->verbose) {
      fprintf(stderr, "%s: hess = \n", me);
      ell3mPRINT(stderr, hess);
    }
  }
  if (1 & (query >> gageSclLaplacian)) {
    ans[offset[gageSclLaplacian]] = hess[0] + hess[4] + hess[8];
    if (ctx->verbose) {
      fprintf(stderr, "%s: lapl = %g + %g + %g  = %g\n", me,
	      hess[0], hess[4], hess[8], ans[offset[gageSclLaplacian]]);
    }
  }
  if (1 & (query >> gageSclHessEval)) {
    ELL_3M_COPY(tmpMat, hess);
    /* HEY: look at the return value for root multiplicity? */
    /* NB: we have solve and then copy because of possible type
       mismatch between double and gage_t */
    ell3mEigensolve(heval, hevec, tmpMat, AIR_TRUE);
    ELL_3V_COPY(ans+offset[gageSclHessEval], heval);
  }
  if (1 & (query >> gageSclHessEvec)) {
    ELL_3M_COPY(ans+offset[gageSclHessEvec], hevec);
  }
  if (1 & (query >> gageScl2ndDD)) {
    ELL_3MV_MUL(tmpVec, hess, norm);
    ans[offset[gageScl2ndDD]] = ELL_3V_DOT(norm, tmpVec);
  }
  if (1 & (query >> gageSclGeomTens)) {
    if (gmag >= ctx->parm.gradMagCurvMin) {
      /* parm.curvNormalSide applied here to determine the sense of the
	 normal when doing all curvature calculations */
      ELL_3M_SCALE(sHess, -(ctx->parm.curvNormalSide)/gmag, hess);
      /* nPerp = I - outer(norm, norm): matrix which projects onto the
	 plane perpendicular to the normal */
      ELL_3MV_OUTER(nProj, norm, norm);
      ELL_3M_SCALE(nPerp, -1, nProj);
      nPerp[0] += 1;
      nPerp[4] += 1;
      nPerp[8] += 1;
      
      /* san->gten = nPerp * sHess * nPerp */
      ELL_3M_MUL(tmpMat, sHess, nPerp);
      ELL_3M_MUL(gten, nPerp, tmpMat);

      if (ctx->verbose) {
	fprintf(stderr, "gten: \n");
	ell3mPRINT(stderr, gten);
	ELL_3MV_MUL(tmpVec, gten, norm);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should be small: %30.15f\n", (double)len);
	ell3vPERP(gp1, norm);
	ELL_3MV_MUL(tmpVec, gten, gp1);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should be bigger: %30.15f\n", (double)len);
	ELL_3V_CROSS(gp2, gp1, norm);
	ELL_3MV_MUL(tmpVec, gten, gp2);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should (also) be bigger: %30.15f\n", (double)len);
      }
    } else {
      ELL_3M_ZERO_SET(gten);
    }
  }
  if (1 && (query >> gageSclCurvedness)) {
    curv = ans[offset[gageSclCurvedness]] = ELL_3M_L2NORM(gten);
  }
  if (1 && (query >> gageSclShapeTrace)) {
    ans[offset[gageSclShapeTrace]] = (curv
				      ? ELL_3M_TRACE(gten)/curv
				      : 0);
  }
  if ( (1 && (query >> gageSclK1)) ||
       (1 && (query >> gageSclK2)) ){
    T = ELL_3M_TRACE(gten);
    N = curv;
    D = 2*N*N - T*T;
    if (D < -0.0000001) {
      fprintf(stderr, "%s: %g %g\n", me, T, N);
      fprintf(stderr, "%s: !!! D curv determinant % 22.10f < 0.0\n", me, D);
      fprintf(stderr, "gten: \n");
      ell3mPRINT(stderr, gten);
    }
    D = AIR_MAX(D, 0);
    D = sqrt(D);
    k1[0] = 0.5*(T + D);
    k2[0] = 0.5*(T - D);
  }
  if (1 && (query >> gageSclMeanCurv)) {
    ans[offset[gageSclMeanCurv]] = (*k1 + *k2)/2;
  }
  if (1 && (query >> gageSclGaussCurv)) {
    ans[offset[gageSclGaussCurv]] = (*k1)*(*k2);
  }
  if (1 && (query >> gageSclShapeIndex)) {
    ans[offset[gageSclShapeIndex]] = -(2/M_PI)*atan2(*k1 + *k2, *k1 - *k2);
  }
  if (1 & (query >> gageSclCurvDir)) {
    /* HEY: this only works when K1, K2, 0 are all well mutually distinct,
       since these are the eigenvalues of the geometry tensor, and this
       code assumes that the eigenspaces are all one-dimensional */
    ELL_3M_COPY(tmpMat, gten);
    ELL_3M_DIAG_SET(tmpMat, gten[0] - *k1, gten[4]- *k1, gten[8] - *k1);
    ell3mNullspace1(tmpVec, tmpMat);
    ELL_3V_COPY(ans+offset[gageSclCurvDir]+0, tmpVec);
    ELL_3M_DIAG_SET(tmpMat, gten[0] - *k2, gten[4] - *k2, gten[8] - *k2);
    ell3mNullspace1(tmpVec, tmpMat);<commit_msg>commented something in or out<commit_after>/*
  teem: <NAME>'s research software
  Copyright (C) 2002, 2001, 2000, 1999, 1998 University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "gage.h"
#include "privateGage.h"

void
_gageSclAnswer (gageContext *ctx, gagePerVolume *pvl) {
  char me[]="_gageSclAnswer";
  unsigned int query;
  gage_t *ans, gmag=0, *hess, *norm, *gvec, *gten, *k1, *k2, sHess[9], curv=0;
  double tmpMat[9], tmpVec[3], hevec[9], heval[3];
  int *offset;

  gage_t len, gp1[3], gp2[3], nPerp[9], nProj[9], ncTen[9];
  double T, N, D;

  query = pvl->query;
  ans = pvl->ans;
  /* convenience pointers for work below */
  offset = gageKindScl->ansOffset;
  hess = ans + offset[gageSclHessian];
  gvec = ans + offset[gageSclGradVec];
  norm = ans + offset[gageSclNormal];
  gten = ans + offset[gageSclGeomTens];
  k1 = ans + offset[gageSclK1];
  k2 = ans + offset[gageSclK2];
  
  if (1 & (query >> gageSclValue)) {
    /* done if doV */
    if (ctx->verbose) {
      fprintf(stderr, "val = % 15.7f\n", (double)(ans[offset[gageSclValue]]));
    }
  }
  if (1 & (query >> gageSclGradVec)) {
    /* done if doD1 */
    if (ctx->verbose) {
      fprintf(stderr, "gvec = ");
      ell3vPRINT(stderr, gvec);
    }
  }
  if (1 & (query >> gageSclGradMag)) {
    /* this is the true value of gradient magnitude */
    gmag = ans[offset[gageSclGradMag]] = sqrt(ELL_3V_DOT(gvec, gvec));
  }

  /* NB: it would seem that gageParmGradMagMin is completely ignored ... */

  if (1 & (query >> gageSclNormal)) {
    if (gmag) {
      ELL_3V_SCALE(norm, 1.0/gmag, gvec);
      /* polishing ... 
      len = sqrt(ELL_3V_DOT(norm, norm));
      ELL_3V_SCALE(norm, 1.0/len, norm);
      */
    } else {
      ELL_3V_COPY(norm, gageZeroNormal);
    }
  }
  if (1 & (query >> gageSclHessian)) {
    /* done if doD2 */
    if (ctx->verbose) {
      fprintf(stderr, "%s: hess = \n", me);
      ell3mPRINT(stderr, hess);
    }
  }
  if (1 & (query >> gageSclLaplacian)) {
    ans[offset[gageSclLaplacian]] = hess[0] + hess[4] + hess[8];
    if (ctx->verbose) {
      fprintf(stderr, "%s: lapl = %g + %g + %g  = %g\n", me,
	      hess[0], hess[4], hess[8], ans[offset[gageSclLaplacian]]);
    }
  }
  if (1 & (query >> gageSclHessEval)) {
    ELL_3M_COPY(tmpMat, hess);
    /* HEY: look at the return value for root multiplicity? */
    /* NB: we have solve and then copy because of possible type
       mismatch between double and gage_t */
    ell3mEigensolve(heval, hevec, tmpMat, AIR_TRUE);
    ELL_3V_COPY(ans+offset[gageSclHessEval], heval);
  }
  if (1 & (query >> gageSclHessEvec)) {
    ELL_3M_COPY(ans+offset[gageSclHessEvec], hevec);
  }
  if (1 & (query >> gageScl2ndDD)) {
    ELL_3MV_MUL(tmpVec, hess, norm);
    ans[offset[gageScl2ndDD]] = ELL_3V_DOT(norm, tmpVec);
  }
  if (1 & (query >> gageSclGeomTens)) {
    if (gmag >= ctx->parm.gradMagCurvMin) {
      /* parm.curvNormalSide applied here to determine the sense of the
	 normal when doing all curvature calculations */
      ELL_3M_SCALE(sHess, -(ctx->parm.curvNormalSide)/gmag, hess);
      /* nPerp = I - outer(norm, norm): matrix which projects onto the
	 plane perpendicular to the normal */
      ELL_3MV_OUTER(nProj, norm, norm);
      ELL_3M_SCALE(nPerp, -1, nProj);
      nPerp[0] += 1;
      nPerp[4] += 1;
      nPerp[8] += 1;
      
      /* san->gten = nPerp * sHess * nPerp */
      ELL_3M_MUL(tmpMat, sHess, nPerp);
      ELL_3M_MUL(gten, nPerp, tmpMat);

      if (ctx->verbose) {
	fprintf(stderr, "gten: \n");
	ell3mPRINT(stderr, gten);
	ELL_3MV_MUL(tmpVec, gten, norm);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should be small: %30.15f\n", (double)len);
	ell3vPERP(gp1, norm);
	ELL_3MV_MUL(tmpVec, gten, gp1);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should be bigger: %30.15f\n", (double)len);
	ELL_3V_CROSS(gp2, gp1, norm);
	ELL_3MV_MUL(tmpVec, gten, gp2);
	len = ELL_3V_LEN(tmpVec);
	fprintf(stderr, "should (also) be bigger: %30.15f\n", (double)len);
      }
    } else {
      ELL_3M_ZERO_SET(gten);
    }
  }
  if (1 && (query >> gageSclCurvedness)) {
    curv = ans[offset[gageSclCurvedness]] = ELL_3M_L2NORM(gten);
  }
  if (1 && (query >> gageSclShapeTrace)) {
    ans[offset[gageSclShapeTrace]] = (curv
				      ? ELL_3M_TRACE(gten)/curv
				      : 0);
  }
  if ( (1 && (query >> gageSclK1)) ||
       (1 && (query >> gageSclK2)) ){
    T = ELL_3M_TRACE(gten);
    N = curv;
    D = 2*N*N - T*T;
    /*
    if (D < -0.0000001) {
      fprintf(stderr, "%s: %g %g\n", me, T, N);
      fprintf(stderr, "%s: !!! D curv determinant % 22.10f < 0.0\n", me, D);
      fprintf(stderr, "gten: \n");
      ell3mPRINT(stderr, gten);
    }
    */
    D = AIR_MAX(D, 0);
    D = sqrt(D);
    k1[0] = 0.5*(T + D);
    k2[0] = 0.5*(T - D);
  }
  if (1 && (query >> gageSclMeanCurv)) {
    ans[offset[gageSclMeanCurv]] = (*k1 + *k2)/2;
  }
  if (1 && (query >> gageSclGaussCurv)) {
    ans[offset[gageSclGaussCurv]] = (*k1)*(*k2);
  }
  if (1 && (query >> gageSclShapeIndex)) {
    ans[offset[gageSclShapeIndex]] = -(2/M_PI)*atan2(*k1 + *k2, *k1 - *k2);
  }
  if (1 & (query >> gageSclCurvDir)) {
    /* HEY: this only works when K1, K2, 0 are all well mutually distinct,
       since these are the eigenvalues of the geometry tensor, and this
       code assumes that the eigenspaces are all one-dimensional */
    ELL_3M_COPY(tmpMat, gten);
    ELL_3M_DIAG_SET(tmpMat, gten[0] - *k1, gten[4]- *k1, gten[8] - *k1);
    ell3mNullspace1(tmpVec, tmpMat);
    ELL_3V_COPY(ans+offset[gageSclCurvDir]+0, tmpVec);
    ELL_3M_DIAG_SET(tmpMat, gten[0] - *k2, gten[4] - *k2, gten[8] - *k2);
    ell3mNullspace1(tmpVec, tmpMat);
{'fraction_non_alphanumeric': 0.15476379831398124, 'fraction_numerical': 0.03340225862891681, 'mean_word_length': 3.137874300756828, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6160460', 'n_tokens_mistral': 5312, 'n_tokens_neox': 5103, 'n_words': 1239}