text:
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <mpfr.h>
#include <petsc.h>
#include "ellipsoid.h"
#include "../sphere/sphere.h"
#include "ellSolv.h"

#undef __FUNCT__
#define __FUNCT__ "InitEllipsoidalAndConvertPoints"
PetscErrorCode InitEllipsoidalAndConvertPoints(EllipsoidalSystem *e, PetscReal a, PetscReal b, PetscReal c, Vec xyz, Vec ell)
{
  PetscErrorCode ierr;
  PetscLogEvent initEvent;
  PetscFunctionBegin;

  ierr = PetscLogEventRegister("Init ellipsoidal system", 0, &initEvent);CHKERRQ(ierr);
  ierr = PetscLogEventBegin(initEvent, 0, 0, 0, 0);CHKERRQ(ierr);
  ierr = initEllipsoidalSystem(e, a, b, c, 32);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(initEvent, 0, 0, 0, 0);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "HowMany"
PetscErrorCode HowMany(PetscInt N, PetscInt *num)
{
  PetscErrorCode ierr;
  PetscInt n, p;
  PetscFunctionBegin;
  *num = 0;
  for(n = 0; n <= N; ++n) {
    for(p=0; p < 2*n + 1; ++p)
      *num += 1;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcEllipsoidFreeEnergy"
PetscErrorCode CalcEllipsoidFreeEnergy(EllipsoidalSystem *e, PetscReal eps1, PetscReal eps2, PetscInt nSrc, Vec srcXYZ, Vec srcMag, PetscReal tol, PetscInt Nmax, Vec tarSol, PetscReal *freeE)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  //EllipsoidalSystem e;
  PetscReal x, y, z;
  PetscInt n, p, npsize;
  Vec srcEll, tarEll;
  PetscScalar *tarEllArray;
  Vec coulCoefs, reactCoefs, extCoefs;
  const PetscScalar *coulCoefsArray, *reactCoefsArray, *extCoefsArray;
  const PetscScalar *srcMagArray;
  Vec EnpVals, FnpVals;
  const PetscScalar *EnpValsArray, *FnpValsArray;
  PetscScalar *tarSolArray;
  PetscFunctionBegin;
  flopCount = 0;

  /* init ellipsoidal system */
  //ierr = initEllipsoidalSystem(&e, a, b, c);CHKERRQ(ierr);

  /* initialize expansion vectors */
  ierr = HowMany(Nmax, &npsize);CHKERRQ(ierr);
  //ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &coulCoefs);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &reactCoefs);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &extCoefs);CHKERRQ(ierr);

  /* initialize point vectors */
  ierr = VecDuplicate(srcXYZ, &srcEll);CHKERRQ(ierr);
  ierr = VecDuplicate(srcEll, &tarEll);CHKERRQ(ierr);

  /* convert points from cartesian to ellipsoidal */
  ierr = CartesianToEllipsoidalVec(e, srcXYZ, srcEll);CHKERRQ(ierr);
  ierr = VecCopy(srcEll, tarEll);CHKERRQ(ierr); //source+target same

  /* Calculate Gnp */
  ierr = CalcCoulombEllCoefs(e, nSrc, srcEll, srcMag, Nmax, &coulCoefs);CHKERRQ(ierr);
  /* calculate Bnp and Cnp */
  ierr = CalcReactAndExtCoefsFromCoulomb(e, eps1, eps2, Nmax, coulCoefs, reactCoefs, extCoefs);

  // create EnpVals and FnpVals vectors
  ierr = VecCreateSeq(PETSC_COMM_SELF, nSrc, &EnpVals);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nSrc, &FnpVals);CHKERRQ(ierr);

  // get read-only pointers for expansion coefficients
  ierr = VecGetArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(reactCoefs, &reactCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(extCoefs, &extCoefsArray);CHKERRQ(ierr);
  // get write pointer for solution vector
  ierr = VecZeroEntries(tarSol);CHKERRQ(ierr);
  ierr = VecGetArray(tarSol, &tarSolArray);CHKERRQ(ierr);
  PetscInt ind = 0;
  for(n=0; n <= Nmax; ++n) {
    for(p=0; p < 2*n+1; ++p) {
      PetscReal Gnp =  coulCoefsArray[ind];
      PetscReal Cnp =   extCoefsArray[ind];
      PetscReal Bnp = reactCoefsArray[ind];

      ierr = CalcSolidInteriorHarmonicVec(e, tarEll, n, p, EnpVals);CHKERRQ(ierr);
      //ierr = CalcSolidExteriorHarmonicVec(&e, tarEll, n, p, FnpVals);CHKERRQ(ierr);

      ierr = VecGetArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      //ierr = VecGetArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);
      for(PetscInt k=0; k < nSrc; ++k) {
	PetscReal Enp = EnpValsArray[k];
	//PetscReal Fnp = FnpValsArray[k];
	tarSolArray[k] += Bnp*Enp; flopCount += 2;
      }
      ierr = VecRestoreArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      //ierr = VecRestoreArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);

      ind++;
    }
  }
  // get read-only pointers for expansion coefficients
  ierr = VecRestoreArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(reactCoefs, &reactCoefsArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(extCoefs, &extCoefsArray);CHKERRQ(ierr);
  *freeE = 0;
  ierr = VecGetArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  /* calculate free energy */
  for(PetscInt k=0; k < nSrc; ++k) {
    printf("sum[d] = %15.15f\n", srcMagArray[k]*tarSolArray[k]);
    *freeE += srcMagArray[k]*tarSolArray[k];
  }
  ierr = VecRestoreArrayRead(srcMag, &srcMagArray);CHKERRQ(ierr);
  ierr = VecRestoreArray(tarSol, &tarSolArray);CHKERRQ(ierr);
  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcEllipsoidTester"
PetscErrorCode CalcEllipsoidTester(PetscReal a, PetscReal b, PetscReal c, PetscReal eps1, PetscReal eps2, PetscInt nSrc, Vec srcXYZ, Vec srcMag, PetscInt nTar, Vec tarXYZ, PetscInt Nmax, Vec tarSol)
{
  PetscErrorCode ierr;
  EllipsoidalSystem e;
  PetscReal x, y, z;
  PetscInt n, p, npsize;
  Vec srcEll, tarEll;
  PetscScalar *tarEllArray;
  Vec coulCoefs, reactCoefs, extCoefs;
  const PetscScalar *coulCoefsArray, *reactCoefsArray, *extCoefsArray;
  Vec EnpVals, FnpVals;
  const PetscScalar *EnpValsArray, *FnpValsArray;
  PetscScalar *tarSolArray;
  PetscFunctionBegin;

  /* init ellipsoidal system */
  ierr = initEllipsoidalSystem(&e, a, b, c, 32);CHKERRQ(ierr);

  /* initialize expansion vectors */
  ierr = HowMany(Nmax, &npsize);CHKERRQ(ierr);
  //ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &coulCoefs);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &reactCoefs);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, npsize, &extCoefs);CHKERRQ(ierr);

  /* initialize ellipsoid vectors */
  ierr = VecDuplicate(srcXYZ, &srcEll);CHKERRQ(ierr);
  ierr = VecDuplicate(tarXYZ, &tarEll);CHKERRQ(ierr);

  /* convert points from cartesian to ellipsoidal */
  ierr = CartesianToEllipsoidalVec(&e, srcXYZ, srcEll);CHKERRQ(ierr);
  ierr = CartesianToEllipsoidalVec(&e, tarXYZ, tarEll);CHKERRQ(ierr);

  /* Calculate Gnp */
  ierr = CalcCoulombEllCoefs(&e, nSrc, srcEll, srcMag, Nmax, &coulCoefs);CHKERRQ(ierr);
  //printf("\n\n##################################\n########## COUL COEFS #############\n###########################################\n\n");
  //printf("nmax: %d\n", Nmax);
  //ierr = VecView(coulCoefs, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  //printf("\n\n##################################\n########## OVER #############\n###########################################\n\n");
  /* calculate Bnp and Cnp */
  ierr = CalcReactAndExtCoefsFromCoulomb(&e, eps1, eps2, Nmax, coulCoefs, reactCoefs, extCoefs);

  // create EnpVals and FnpVals vectors
  ierr = VecCreateSeq(PETSC_COMM_SELF, nTar, &EnpVals);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nTar, &FnpVals);CHKERRQ(ierr);

  // get read-only pointers for expansion coefficients
  ierr = VecGetArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(reactCoefs, &reactCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(extCoefs, &extCoefsArray);CHKERRQ(ierr);
  // get write pointer for solution vector
  ierr = VecZeroEntries(tarSol);CHKERRQ(ierr);
  ierr = VecGetArray(tarSol, &tarSolArray);CHKERRQ(ierr);
  PetscInt ind = 0;
  for(n=0; n <= Nmax; ++n) {
    printf("n: %d\n", n);
    for(p=0; p < 2*n+1; ++p) {
      PetscReal Gnp =  coulCoefsArray[ind];
      PetscReal Cnp =   extCoefsArray[ind];
      PetscReal Bnp = reactCoefsArray[ind];

      ierr = CalcSolidInteriorHarmonicVec(&e, tarEll, n, p, EnpVals);CHKERRQ(ierr);
      ierr = CalcSolidExteriorHarmonicVec(&e, tarEll, n, p, FnpVals);CHKERRQ(ierr);

      ierr = VecGetArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      ierr = VecGetArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);
      for(PetscInt k=0; k < nTar; ++k) {
	if(n==0 && p==0) {
	  tarSolArray[k] = 0;
	}
	PetscReal lambda;
	PetscInt index = 3*k+0;
	ierr = VecGetValues(tarEll, 1, &index, &lambda);CHKERRQ(ierr);

	PetscReal Enp = EnpValsArray[k];
	PetscReal Fnp = FnpValsArray[k];
	if(PetscAbsReal(lambda) <= a)
	  tarSolArray[k] += Bnp*Enp; //(Gnp/eps1)*Fnp;
	else {
	  tarSolArray[k] += Cnp*Fnp;//(Gnp/eps1)*Fnp;
	}

      }
      ierr = VecRestoreArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      ierr = VecRestoreArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);

      ind++;
    }
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcEllipsoidCoulombPotential"
PetscErrorCode CalcEllipsoidCoulombPotential(PetscReal a, PetscReal b, PetscReal c, PetscReal eps1, PetscReal eps2, PetscInt nCharges, Vec chargeXYZ, Vec chargeMag, PetscInt nSol, Vec solXYZ, PetscInt Nmax, Vec targetSol)
{
  PetscErrorCode ierr;
  Vec chargeEll;
  Vec solEll;
  Vec coulCoefs;
  Vec FnpVals;
  EllipsoidalSystem e;
  PetscScalar *vecPtr;
  PetscScalar *coulCoefsArray;
  PetscScalar *FnpValsArray;
  PetscFunctionBegin;

  ierr = initEllipsoidalSystem(&e, a, b, c, 32);CHKERRQ(ierr);

  // create charge ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nCharges, &chargeEll);CHKERRQ(ierr);
  printf("Converting charge points to ellipsoidal\n");
  ierr = CartesianToEllipsoidalVec(&e, chargeXYZ, chargeEll);CHKERRQ(ierr);
  // create solution ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nSol, &solEll);CHKERRQ(ierr);
  printf("Converting solution points to ellipsoidal\n");
  ierr = CartesianToEllipsoidalVec(&e, solXYZ, solEll);CHKERRQ(ierr);

  //calculate coulomb coefficients
  ierr = CalcCoulombEllCoefs(&e, nCharges, chargeEll, chargeMag, Nmax, &coulCoefs);CHKERRQ(ierr);

  //initialize vector for interior harmonic calculations
  ierr = VecCreateSeq(PETSC_COMM_SELF, nSol, &FnpVals);CHKERRQ(ierr);

  ierr = VecGetArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArray(targetSol, &vecPtr);CHKERRQ(ierr);
  PetscInt index = 0;
  for(PetscInt n=0; n<=Nmax; ++n) {
    for(PetscInt p=0; p < 2*n+1; ++p) {
      //calculate Enp Vals
      ierr = CalcSolidExteriorHarmonicVec(&e, solEll, n, p, FnpVals);CHKERRQ(ierr);
      ierr = VecGetArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);
      for(PetscInt k=0; k<nSol; ++k) {
	vecPtr[k] += coulCoefsArray[index]*FnpValsArray[k];
      }
      ierr = VecRestoreArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);
      index++;
    }
  }
  ierr = VecRestoreArray(targetSol, &vecPtr);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);

  ierr = VecDestroy(&solEll);CHKERRQ(ierr);
  ierr = VecDestroy(&chargeEll);CHKERRQ(ierr);
  ierr = VecDestroy(&coulCoefs);CHKERRQ(ierr);
  ierr = VecDestroy(&FnpVals);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcEllipsoidSolvationPotential"
PetscErrorCode CalcEllipsoidSolvationPotential(PetscReal a, PetscReal b, PetscReal c, PetscReal eps1, PetscReal eps2, PetscInt nSource, Vec sourceXYZ, Vec sourceMag, PetscInt nTarget, Vec targetXYZ, PetscInt Nmax, Vec targetSol)
{
  PetscErrorCode ierr;
  EllipsoidalSystem e;
  PetscReal x, y, z;
  PetscInt n, p;
  PetscInt intPts, extPts;
  Vec tarIntXYZ, tarExtXYZ;
  PetscScalar *tarIntXYZArray, *tarExtXYZArray;
  Vec tarIntEll, tarExtEll;
  Vec srcEll, tarEll;
  Vec coulCoefs, reactCoefs, extCoefs;
  Vec EnpVals, FnpVals;
  PetscInt *isExt;
  const PetscScalar *targetXYZArray, *coulCoefsArray, *reactCoefsArray, *extCoefsArray;
  PetscScalar *targetSolArray;
  const PetscScalar *EnpValsArray, *FnpValsArray;

  //PetscLogStage stageIntSolids, stageExtSolids;
  PetscLogEvent *INTERIOR_FLOPS;
  PetscLogEvent *EXTERIOR_FLOPS;
  //PetscLogDouble user_log_flops;
  PetscFunctionBegin;

  ierr = initEllipsoidalSystem(&e, a, b, c, 32);CHKERRQ(ierr);

  // calculate the number of interior and exterior points
  ierr = VecGetArrayRead(targetXYZ, &targetXYZArray);CHKERRQ(ierr);
  intPts = 0; extPts = 0;
  for(int k=0; k < nTarget; ++k) {
    x = targetXYZArray[3*k+0];
    y = targetXYZArray[3*k+1];
    z = targetXYZArray[3*k+2];
    if((x*x)/(a*a) + (y*y)/(b*b) + (z*z)/(c*c) <= 1)
      intPts++;
    else
      extPts++;
  }
  // init interior and exterior target point vectors
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*intPts, &tarIntXYZ);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*extPts, &tarExtXYZ);CHKERRQ(ierr);
  // sort points into int and ext vectors
  ierr = VecGetArray(tarIntXYZ, &tarIntXYZArray);CHKERRQ(ierr);
  ierr = VecGetArray(tarExtXYZ, &tarExtXYZArray);CHKERRQ(ierr);

  ierr = PetscMalloc1(sizeof(PetscInt)*nTarget, &isExt);CHKERRQ(ierr);

  intPts = 0; extPts = 0;
  for(int k=0; k < nTarget; ++k) {
    x = targetXYZArray[3*k+0];
    y = targetXYZArray[3*k+1];
    z = targetXYZArray[3*k+2];
    if((x*x)/(a*a) + (y*y)/(b*b) + (z*z)/(c*c) <= 1) {
      tarIntXYZArray[3*intPts+0] = x;
      tarIntXYZArray[3*intPts+1] = y;
      tarIntXYZArray[3*intPts+2] = z;
      isExt[k] = 0;
      intPts++;
    }
    else {
      tarExtXYZArray[3*extPts+0] = x;
      tarExtXYZArray[3*extPts+1] = y;
      tarExtXYZArray[3*extPts+2] = z;
      isExt[k] = 1;
      extPts++;
    }
  }
  ierr = VecRestoreArray(tarIntXYZ, &tarIntXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArray(tarExtXYZ, &tarExtXYZArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(targetXYZ, &targetXYZArray);CHKERRQ(ierr);

  // create source ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nSource, &srcEll);CHKERRQ(ierr);
  ierr = CartesianToEllipsoidalVec(&e, sourceXYZ, srcEll);CHKERRQ(ierr);
  // create target ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nTarget, &tarEll);CHKERRQ(ierr);
  ierr = CartesianToEllipsoidalVec(&e, targetXYZ, tarEll);CHKERRQ(ierr);
  // create target interior ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*intPts, &tarIntEll);CHKERRQ(ierr);
  ierr = CartesianToEllipsoidalVec(&e, tarIntXYZ, tarIntEll);CHKERRQ(ierr);
  // create target exterior ellipsoidal vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*extPts, &tarExtEll);CHKERRQ(ierr);
  ierr = CartesianToEllipsoidalVec(&e, tarExtXYZ, tarExtEll);CHKERRQ(ierr);

  // calculate coulomb coefficients
  ierr = CalcCoulombEllCoefs(&e, nSource, srcEll, sourceMag, Nmax, &coulCoefs);CHKERRQ(ierr);
  // init reacCoefs, extCoefs from size of coulCoefs
  ierr = VecDuplicate(coulCoefs, &reactCoefs);CHKERRQ(ierr);
  ierr = VecDuplicate(coulCoefs, &extCoefs);CHKERRQ(ierr);
  // calc reaction and exterior expansion coefs from coulomb coefs
  ierr = CalcReactAndExtCoefsFromCoulomb(&e, eps1, eps2, Nmax, coulCoefs, reactCoefs, extCoefs);

  // create EnpVals and FnpVals vectors
  ierr = VecCreateSeq(PETSC_COMM_SELF, intPts, &EnpVals);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, extPts, &FnpVals);CHKERRQ(ierr);

  // get read-only pointers for expansion coefficients
  ierr = VecGetArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(reactCoefs, &reactCoefsArray);CHKERRQ(ierr);
  ierr = VecGetArrayRead(extCoefs, &extCoefsArray);CHKERRQ(ierr);
  // get write pointer for solution vector
  ierr = VecZeroEntries(targetSol);CHKERRQ(ierr);
  ierr = VecGetArray(targetSol, &targetSolArray);CHKERRQ(ierr);
  // loop
  PetscInt howmany = 0;
  for(n=0; n <= Nmax; ++n) {
    for(p=0; p < 2*n+1; ++p)
      howmany++;
  }

  ierr = PetscMalloc1(sizeof(PetscLogEvent)*howmany, &INTERIOR_FLOPS);CHKERRQ(ierr);
  ierr = PetscMalloc1(sizeof(PetscLogEvent)*howmany, &EXTERIOR_FLOPS);CHKERRQ(ierr);
  char intStr[40] = "Interior Harmonic %d";
  char extStr[40] = "Exterior Harmonic %d";
  char indvIntStr[40];
  char indvExtStr[40];
  for(n=0; n<howmany; ++n) {
    sprintf(indvIntStr, intStr, n);
    sprintf(indvExtStr, extStr, n);
    ierr = PetscLogEventRegister(indvIntStr, 0, INTERIOR_FLOPS+n);CHKERRQ(ierr);
    ierr = PetscLogEventRegister(indvExtStr, 0, EXTERIOR_FLOPS+n);CHKERRQ(ierr);
  }
  PetscInt ind = 0;
  PetscInt intInd;
  PetscInt extInd;
  for(n=0; n <= Nmax; ++n) {
    printf("%d/%d\n", n, Nmax);
    for(p=0; p < 2*n+1; ++p) {

      PetscReal Gnp =  coulCoefsArray[ind];
      PetscReal Cnp =   extCoefsArray[ind];
      PetscReal Bnp = reactCoefsArray[ind];
      // calc solid harmonics for interior,exterior points

      //ierr = PetscLogStagePush(stageIntSolids);CHKERRQ(ierr);
      ierr = PetscLogEventBegin(INTERIOR_FLOPS[ind],0,0,0,0);CHKERRQ(ierr);
      ierr = CalcSolidInteriorHarmonicVec(&e, tarIntEll, n, p, EnpVals);CHKERRQ(ierr);
      ierr = PetscLogEventEnd  (INTERIOR_FLOPS[ind],0,0,0,0);CHKERRQ(ierr);
      //ierr = PetscLogStagePop();

      //ierr = PetscLogStagePush(stageExtSolids);CHKERRQ(ierr);
      ierr = PetscLogEventBegin(EXTERIOR_FLOPS[ind],0,0,0,0);CHKERRQ(ierr);
      ierr = CalcSolidExteriorHarmonicVec(&e, tarExtEll, n, p, FnpVals);CHKERRQ(ierr);
      ierr = PetscLogEventEnd(EXTERIOR_FLOPS[ind],0,0,0,0);CHKERRQ(ierr);
      //ierr = PetscLogStagePop();CHKERRQ(ierr);

      ierr = VecGetArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      ierr = VecGetArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);
      intInd = 0;
      extInd = 0;
      for(PetscInt k=0; k<nTarget; ++k) {

	if(isExt[k] == 0) { // interior points
	  //printf("E%d%d: %15.15f\nBnp %d%d: %15.15f\n\n", n, p, n, p, EnpValsArray[intInd], reactCoefsArray[ind]);
	  targetSolArray[k] += reactCoefsArray[ind]*EnpValsArray[intInd];
	  //ierr = CalcSolidInterior
	  //targetSolArray[k] += Gnp*EnpValsArray[intInd];
	  //targetSolArray[k] += 1;
	  intInd++;
	}
	else { // exterior points
	  //printf("F%d%d: %15.15f\nCnp %d%d: %15.15f\n\n", n, p, n, p, FnpValsArray[intInd], extCoefsArray[ind]);
	  //targetSolArray[k] += extCoefsArray[ind]*FnpValsArray[extInd];
	  //targetSolArray[k] += (Gnp*FnpValsArray[extInd])/eps1;
	  targetSolArray[k] += 0;
	  extInd++;
	}
      }
      ierr = VecRestoreArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      ierr = VecRestoreArrayRead(FnpVals, &FnpValsArray);CHKERRQ(ierr);

      ind++;

    }
  }
  //ierr = PetscLogEventEnd(SOLID_FLOPS, 0, 0, 0, 0);

  ierr = VecRestoreArray(targetSol, &targetSolArray);CHKERRQ(ierr);
  // restore read-only pointers for expansion coefficients
  ierr = VecRestoreArrayRead(coulCoefs , &coulCoefsArray );CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(reactCoefs, &reactCoefsArray);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(extCoefs  , &extCoefsArray  );CHKERRQ(ierr);

  //output flops to file
  FILE *fpext = fopen("flopsext.txt", "w");
  FILE *fpint = fopen("flopsint.txt", "w");
  PetscStageLog stageLog;
  ierr = PetscLogGetStageLog(&stageLog);CHKERRQ(ierr);

  /*
  for(int i=0; i<howmany; ++i) {
    fprintf(fpext, "Ext %d = %.4e\n", i, stageLog->stageInfo[stageExtSolids].eventLog->eventInfo[EXTERIOR_FLOPS[i]].flops);
    fprintf(fpint, "Int %d = %.4e\n", i, stageLog->stageInfo[stageIntSolids].eventLog->eventInfo[INTERIOR_FLOPS[i]].flops);
  }
  */

  ierr = VecDestroy(&tarIntXYZ);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&tarExtXYZ);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&tarIntEll);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&tarExtEll);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&srcEll);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&tarEll);CHKERRQ(ierr);CHKERRQ(ierr);
  ierr = VecDestroy(&coulCoefs);CHKERRQ(ierr);
  ierr = VecDestroy(&reactCoefs);CHKERRQ(ierr);
  ierr = VecDestroy(&extCoefs);CHKERRQ(ierr);
  ierr = VecDestroy(&EnpVals); CHKERRQ(ierr);
  ierr = VecDestroy(&FnpVals); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcSolidInteriorHarmonic"
PetscErrorCode CalcSolidInteriorHarmonic(EllipsoidalSystem* e, PetscReal lambda, PetscReal mu, PetscReal nu, PetscInt n, PetscInt p, PetscReal *val)
{
  PetscErrorCode ierr;
  PetscInt signm = 1;
  PetscInt signn = 1;
  PetscReal eL, eM, eN;
  PetscFunctionBegin;

  if(mu < 0)
    signm = -1;
  if(nu < 0)
    signn = -1;

  calcLame(e, n, p, lambda, signm, signn, &eL);
  calcLame(e, n, p, mu, signm, signn, &eM);
  calcLame(e, n, p, nu, signm, signn, &eN);

  *val = eL*eM*eN;

  ierr= PetscLogFlops(3);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcSolidInteriorHarmonicVec"
/*
  ellPoints is ellPoints of size 3*nPoints containing ellipsoidal coordinates
  outputs solid harmonics to values of size nPoints
*/
PetscErrorCode CalcSolidInteriorHarmonicVec(EllipsoidalSystem* e, Vec ellPoints, PetscInt n, PetscInt p, Vec values)
{
  PetscErrorCode ierr;
  PetscInt signm = 1;
  PetscInt signn = 1;
  PetscInt nPoints;
  PetscReal lambda, mu, nu;
  PetscReal eL, eM, eN;
  PetscReal Enp;
  const PetscScalar* ellPointsArray;
  PetscFunctionBegin;
  ierr = VecGetSize(ellPoints, &nPoints);CHKERRQ(ierr);
  nPoints = nPoints/3;

  ierr = VecGetArrayRead(ellPoints, &ellPointsArray); CHKERRQ(ierr);
  for(PetscInt k=0; k<nPoints; ++k) {
    lambda = ellPointsArray[3*k+0];
    mu     = ellPointsArray[3*k+1];
    nu     = ellPointsArray[3*k+2];
    if(mu < 0)
      signm = -1;
    else signm = 1;
    if(nu < 0)
      signn = -1;
    else signn = 1;
    ierr = calcLame(e, n, p, lambda, signm, signn, &eL);
    ierr = calcLame(e, n, p, mu, signm, signn, &eM);
    ierr = calcLame(e, n, p, nu, signm, signn, &eN);
    Enp = eL*eM*eN;
    ierr = VecSetValue(values, k, Enp, INSERT_VALUES); CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(values); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(values); CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(ellPoints, &ellPointsArray); CHKERRQ(ierr);

  ierr = PetscLogFlops(3*nPoints);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcSolidExteriorHarmonic"
PetscErrorCode CalcSolidExteriorHarmonic(EllipsoidalSystem *e, PetscReal lambda, PetscReal mu, PetscReal nu, PetscInt n, PetscInt p, PetscReal *val)
{
  PetscErrorCode ierr;
  PetscReal Enp, Inp;
  PetscInt signm;
  PetscInt signn;
  PetscFunctionBegin;

  if(mu < 0)
    signm = -1;
  else signm = 1;
  if(nu < 0)
    signn = -1;
  else signn = 1;

  ierr = CalcSolidInteriorHarmonic(e, lambda, mu, nu, n, p, &Enp);CHKERRQ(ierr);
  ierr = calcI(e, n, p, lambda, signm, signn, &Inp);CHKERRQ(ierr);

  *val = (2*n + 1) * Enp * Inp;

  PetscLogFlops(4);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcSolidExteriorHarmonicVec"
PetscErrorCode CalcSolidExteriorHarmonicVec(EllipsoidalSystem* e, Vec ellPoints, PetscInt n, PetscInt p, Vec values)
{
  PetscErrorCode ierr;
  PetscReal Enp, Inp, Fnp;
  PetscReal lambda, mu, nu;
  PetscInt signm;
  PetscInt signn;
  PetscInt nPoints;
  const PetscScalar* ellPointsArray;
  PetscFunctionBegin;
  ierr = VecGetSize(ellPoints, &nPoints);CHKERRQ(ierr);
  nPoints = nPoints/3;

  ierr = VecGetArrayRead(ellPoints, &ellPointsArray);CHKERRQ(ierr);
  for(PetscInt k=0; k<nPoints; ++k) {
    lambda = ellPointsArray[3*k+0];
    mu = ellPointsArray[3*k+1];
    nu = ellPointsArray[3*k+2];

    if(mu < 0)
      signm = -1;
    else signm = 1;
    if(nu < 0)
      signn = -1;
    else signn = 1;

    ierr = CalcSolidInteriorHarmonic(e, lambda, mu, nu, n, p, &Enp);CHKERRQ(ierr);
    ierr = calcI(e, n, p, lambda, signm, signn, &Inp);CHKERRQ(ierr);
    Fnp = (2*n + 1) * Enp * Inp;
    ierr = VecSetValues(values, 1, &k, &Fnp, INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(values);CHKERRQ(ierr); ierr = VecAssemblyEnd(values);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(ellPoints, &ellPointsArray);CHKERRQ(ierr);

  ierr = PetscLogFlops(4*nPoints);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcCoulombEllCoefs"
/*
  GnpVals should be uninitialized. I
*/
PetscErrorCode CalcCoulombEllCoefs(EllipsoidalSystem* e, PetscInt nSource, Vec srcPoints, Vec srcCharges, PetscInt Nmax, Vec* GnpVals)
{
  PetscErrorCode ierr;
  PetscReal normConstant;
  PetscReal Gnp;
  PetscInt count;

  Vec EnpVals;
  const PetscScalar* EnpValsArray;
  const PetscScalar* srcPointsArray;
  const PetscScalar* srcChargesArray;
  PetscFunctionBegin;

  ierr = VecCreateSeq(PETSC_COMM_SELF, nSource, &EnpVals);CHKERRQ(ierr);

  count = 0;
  for(PetscInt n=0; n<=Nmax; ++n) {
    for(PetscInt p=0; p<2*n+1; ++p)
      count++;
  }
  //ierr = VecDestroy(GnpVals);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, count, GnpVals);CHKERRQ(ierr);

  ierr = VecGetArrayRead(srcCharges, &srcChargesArray);CHKERRQ(ierr);
  count = 0;
  for(PetscInt n=0; n<=Nmax; ++n) {
    for(PetscInt p=0; p<2*n+1; ++p) {
      ierr = calcNormalization(e, n, p, &normConstant);CHKERRQ(ierr);

      ierr = CalcSolidInteriorHarmonicVec(e, srcPoints, n, p, EnpVals);

      ierr = VecGetArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      Gnp = 0;
      for(PetscInt k=0; k<nSource; ++k) {
	Gnp += srcChargesArray[k]*EnpValsArray[k];
      }
      ierr = VecRestoreArrayRead(EnpVals, &EnpValsArray);CHKERRQ(ierr);
      Gnp *= (4*PETSC_PI)/((2.0*n+1.0)*normConstant);
      ierr = VecSetValues(*GnpVals, 1, &count, &Gnp, INSERT_VALUES);CHKERRQ(ierr);
      count++;
    }
  }
  ierr = VecAssemblyBegin(*GnpVals);CHKERRQ(ierr); ierr = VecAssemblyEnd(*GnpVals);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(srcCharges, &srcChargesArray);CHKERRQ(ierr);

  ierr = VecDestroy(&EnpVals); CHKERRQ(ierr);

  ierr = PetscLogFlops(count * (6 + 3*nSource));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "CalcReactAndExtCoefsFromCoulomb"
PetscErrorCode CalcReactAndExtCoefsFromCoulomb(EllipsoidalSystem* e, PetscReal eps1, PetscReal eps2, PetscInt Nmax, Vec coulCoefs, Vec reacCoefs, Vec extCoefs)
{
  PetscErrorCode ierr;
  PetscReal Ea   , Ia   , Fa;
  PetscReal EaDer, IaDer, FaDer;
  PetscReal Gnp, Bnp, Cnp;
  PetscReal temp;
  const PetscScalar* coulCoefsArray;
  PetscInt count;
  PetscFunctionBegin;

  ierr = VecGetArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);
  count = 0;
  for(PetscInt n=0; n<=Nmax; ++n) {
    for(PetscInt p=0; p<2*n+1; ++p) {

      Gnp = coulCoefsArray[count];

      ierr = calcLame(e, n, p, e->a, 1, 1, &Ea);CHKERRQ(ierr);
      ierr = calcI   (e, n, p, e->a, 1, 1, &Ia);CHKERRQ(ierr);
      Fa    = (2*n + 1) * Ea    * Ia;
      calcLameDerivative(e, n, p, e->a, 1, 1, &EaDer);
      calcIDerivative   (e, n, p, e->a, 1, 1, &IaDer);
      FaDer = (2*n + 1) * (Ea*IaDer + EaDer*Ia);

      temp = (Fa/Ea)*(eps1 - eps2)/(eps1*eps2);
      temp /= (1 - (eps1/eps2)*((EaDer*Fa)/(FaDer*Ea)));
      Bnp  = temp*Gnp;

      Cnp = (eps1/eps2)*(EaDer/FaDer)*Bnp;
      Cnp += (Gnp/eps2);

      ierr = VecSetValues(reacCoefs, 1, &count, &Bnp, INSERT_VALUES);CHKERRQ(ierr);
      ierr = VecSetValues(extCoefs, 1, &count, &Cnp, INSERT_VALUES);CHKERRQ(ierr);

      count++;
    }
  }

  ierr = VecAssemblyBegin(reacCoefs);CHKERRQ(ierr);
  ierr = VecAssemblyEnd  (reacCoefs);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(extCoefs);CHKERRQ(ierr);
  ierr = VecAssemblyEnd  (extCoefs);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(coulCoefs, &coulCoefsArray);CHKERRQ(ierr);

  ierr = PetscLogFlops(29*count);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:45:10 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2hc-backward 4 */

/*
 * This function contains 34 FP additions, 18 FP multiplications,
 * (or, 28 additions, 12 multiplications, 6 fused multiply/add),
 * 15 stack variables, and 32 memory accesses
 */
static const fftw_real K1_414213562 = FFTW_KONST(+1.414213562373095048801688724209698078569671875);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2hc_backward_4(fftw_real *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_real *X;
     fftw_real *Y;
     X = A;
     Y = A + (4 * iostride);
     {
	  fftw_real tmp39;
	  fftw_real tmp42;
	  fftw_real tmp37;
	  fftw_real tmp40;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp38;
	       fftw_real tmp41;
	       fftw_real tmp35;
	       fftw_real tmp36;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp38 = X[iostride];
	       tmp39 = K2_000000000 * tmp38;
	       tmp41 = Y[-iostride];
	       tmp42 = K2_000000000 * tmp41;
	       tmp35 = X[0];
	       tmp36 = X[2 * iostride];
	       tmp37 = tmp35 + tmp36;
	       tmp40 = tmp35 - tmp36;
	  }
	  X[2 * iostride] = tmp37 - tmp39;
	  X[0] = tmp37 + tmp39;
	  X[3 * iostride] = tmp40 + tmp42;
	  X[iostride] = tmp40 - tmp42;
     }
     X = X + dist;
     Y = Y - dist;
     for (i = 2; i < m; i = i + 2, X = X + dist, Y = Y - dist, W = W + 3) {
	  fftw_real tmp9;
	  fftw_real tmp28;
	  fftw_real tmp18;
	  fftw_real tmp25;
	  fftw_real tmp12;
	  fftw_real tmp24;
	  fftw_real tmp21;
	  fftw_real tmp29;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp7;
	       fftw_real tmp8;
	       fftw_real tmp16;
	       fftw_real tmp17;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp7 = X[0];
	       tmp8 = Y[-2 * iostride];
	       tmp9 = tmp7 + tmp8;
	       tmp28 = tmp7 - tmp8;
	       tmp16 = Y[0];
	       tmp17 = X[2 * iostride];
	       tmp18 = tmp16 - tmp17;
	       tmp25 = tmp16 + tmp17;
	  }
	  {
	       fftw_real tmp10;
	       fftw_real tmp11;
	       fftw_real tmp19;
	       fftw_real tmp20;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp10 = X[iostride];
	       tmp11 = Y[-3 * iostride];
	       tmp12 = tmp10 + tmp11;
	       tmp24 = tmp10 - tmp11;
	       tmp19 = Y[-iostride];
	       tmp20 = X[3 * iostride];
	       tmp21 = tmp19 - tmp20;
	       tmp29 = tmp19 + tmp20;
	  }
	  X[0] = tmp9 + tmp12;
	  {
	       fftw_real tmp14;
	       fftw_real tmp22;
	       fftw_real tmp13;
	       fftw_real tmp15;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp14 = tmp9 - tmp12;
	       tmp22 = tmp18 - tmp21;
	       tmp13 = c_re(W[1]);
	       tmp15 = c_im(W[1]);
	       X[2 * iostride] = (tmp13 * tmp14) + (tmp15 * tmp22);
	       Y[-iostride] = (tmp13 * tmp22) - (tmp15 * tmp14);
	  }
	  Y[-3 * iostride] = tmp18 + tmp21;
	  {
	       fftw_real tmp26;
	       fftw_real tmp30;
	       fftw_real tmp23;
	       fftw_real tmp27;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp26 = tmp24 + tmp25;
	       tmp30 = tmp28 - tmp29;
	       tmp23 = c_re(W[0]);
	       tmp27 = c_im(W[0]);
	       Y[-2 * iostride] = (tmp23 * tmp26) - (tmp27 * tmp30);
	       X[iostride] = (tmp27 * tmp26) + (tmp23 * tmp30);
	  }
	  {
	       fftw_real tmp32;
	       fftw_real tmp34;
	       fftw_real tmp31;
	       fftw_real tmp33;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp32 = tmp25 - tmp24;
	       tmp34 = tmp28 + tmp29;
	       tmp31 = c_re(W[2]);
	       tmp33 = c_im(W[2]);
	       Y[0] = (tmp31 * tmp32) - (tmp33 * tmp34);
	       X[3 * iostride] = (tmp33 * tmp32) + (tmp31 * tmp34);
	  }
     }
     if (i == m) {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp3;
	  fftw_real tmp4;
	  fftw_real tmp5;
	  fftw_real tmp6;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = X[0];
	  tmp2 = X[iostride];
	  tmp3 = tmp1 - tmp2;
	  tmp4 = Y[0];
	  tmp5 = Y[-iostride];
	  tmp6 = tmp4 + tmp5;
	  X[0] = K2_000000000 * (tmp1 + tmp2);
	  X[2 * iostride] = -(K2_000000000 * (tmp4 - tmp5));
	  X[iostride] = K1_414213562 * (tmp3 - tmp6);
	  X[3 * iostride] = -(K1_414213562 * (tmp3 + tmp6));
     }
}

static const int twiddle_order[] =
{1, 2, 3};
fftw_codelet_desc fftw_hc2hc_backward_4_desc =
{
     "fftw_hc2hc_backward_4",
     (void (*)()) fftw_hc2hc_backward_4,
     4,
     FFTW_BACKWARD,
     FFTW_HC2HC,
     102,
     3,
     twiddle_order,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2016-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json_fwd.hpp>
#include <cstdint>
#include <string>

namespace cb::crypto {
enum class Algorithm { MD5, SHA1, SHA256, SHA512 };

bool isSupported(Algorithm algorithm);

const int MD5_DIGEST_SIZE = 16;
const int SHA1_DIGEST_SIZE = 20;
const int SHA256_DIGEST_SIZE = 32;
const int SHA512_DIGEST_SIZE = 64;

/**
 * Generate a HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string HMAC(Algorithm algorithm,
                 std::string_view key,
                 std::string_view data);

/**
 * Generate a PBKDF2_HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string PBKDF2_HMAC(Algorithm algorithm,
                        std::string_view pass,
                        std::string_view salt,
                        unsigned int iterationCount);

/**
 * Generate a digest by using the requested algorithm
 */
std::string digest(Algorithm algorithm, std::string_view data);

} // namespace cb::crypto

ext: h


text:
/* movstat/movmad.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

/*
gsl_movstat_mad()
  Apply a moving MAD to an input vector (with scale factor to make an
unbiased estimate of sigma)

Inputs: endtype - how to handle end points
        x       - input vector, size n
        xmedian - (output) vector of median values of x, size n
                  xmedian_i = median of window centered on x_i
        xmad    - (output) vector of estimated standard deviations of x, size n
                  xmad_i = MAD of i-th window: 1.4826 * median(|x_i - xmedian_i|)
        w       - workspace
*/

int
gsl_movstat_mad(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian, gsl_vector * xmad,
                gsl_movstat_workspace * w)
{
  if (x->size != xmedian->size)
    {
      GSL_ERROR("x and xmedian vectors must have same length", GSL_EBADLEN);
    }
  else if (x->size != xmad->size)
    {
      GSL_ERROR("x and xmad vectors must have same length", GSL_EBADLEN);
    }
  else
    {
      double scale = 1.482602218505602;
      int status = gsl_movstat_apply_accum(endtype, x, gsl_movstat_accum_mad, &scale, xmedian, xmad, w);
      return status;
    }
}

/*
gsl_movstat_mad0()
  Apply a moving MAD to an input vector (without scale factor to make an
unbiased estimate of sigma)

Inputs: endtype - how to handle end points
        x       - input vector, size n
        xmedian - (output) vector of median values of x, size n
                  xmedian_i = median of window centered on x_i
        xmad    - (output) vector of estimated standard deviations of x, size n
                  xmad_i = MAD of i-th window: median(|x_i - xmedian_i|)
        w       - workspace
*/

int
gsl_movstat_mad0(const gsl_movstat_end_t endtype, const gsl_vector * x, gsl_vector * xmedian, gsl_vector * xmad,
                 gsl_movstat_workspace * w)
{
  if (x->size != xmedian->size)
    {
      GSL_ERROR("x and xmedian vectors must have same length", GSL_EBADLEN);
    }
  else if (x->size != xmad->size)
    {
      GSL_ERROR("x and xmad vectors must have same length", GSL_EBADLEN);
    }
  else
    {
      double scale = 1.0;
      int status = gsl_movstat_apply_accum(endtype, x, gsl_movstat_accum_mad, &scale, xmedian, xmad, w);
      return status;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
//
// Created by root on 14.03.20.
//

#ifndef NETWORK_NETWORK_H
#define NETWORK_NETWORK_H

#include <memory>
#include <vector>
#include <functional>
#include <cblas.h>
#include <iostream>
#include <random>
#include "cmake-build-debug/data.pb.h"
#include "Utils.h"

extern std::mt19937 generator;

struct Layer {
    std::unique_ptr<float[]> weights;
    std::unique_ptr<float[]> biases;
    std::unique_ptr<float[]> pre_activ;
    std::unique_ptr<float[]> post_activ;
    std::unique_ptr<float[]> deltas;
    size_t num_rows;
    size_t num_cols;
    std::function<void()> activation;
    std::function<void(float *, float *)> derivative;

    void forward_pass(float *input);

};

class Network {

private:

    std::vector<Layer> layers;
    size_t input_size;

    //looking over backprop implementation
    //Something is not working correctly
    template<typename Loss>
    void back_propagate(float *act_output) {
        {
            Layer &current = get_output_layer();
            const size_t out_neurons = current.num_rows;

            std::unique_ptr<float[]> loss_deriv = std::make_unique<float[]>(out_neurons);
            //Computing the derivative of the loss function with respect to the net_output
            Loss::apply_deriv(current.post_activ.get(), current.post_activ.get() + out_neurons,
                              act_output, loss_deriv.get());

            current.derivative(current.pre_activ.get(), current.deltas.get());
            std::transform(current.deltas.get(), current.deltas.get() + out_neurons, loss_deriv.get(),
                           current.deltas.get(), std::multiplies<float>{});

        }

        //stuff for all the other layers
        for (int l = layers.size() - 2; l >= 0; --l) {
            Layer &last = layers[l + 1];
            Layer &current = layers[l];

            std::unique_ptr<float[]> deriv_out = std::make_unique<float[]>(current.num_rows);
            //
            current.derivative(current.pre_activ.get(), deriv_out.get());

            //now we can compute the deltas of all the other layers
            cblas_sgemv(CblasRowMajor, CblasTrans, last.num_rows, last.num_cols, 1.0f, last.weights.get(),
                        last.num_cols, last.deltas.get(),
                        1, 0.0f, current.deltas.get(), 1);

            std::transform(current.deltas.get(), current.deltas.get() + current.num_rows, deriv_out.get(),
                           current.deltas.get(), std::multiplies<float>{});

        }
    }

public:

    Network(size_t input) : input_size(input) {}

    template<typename Activ>
    void add_layer(size_t num_hidden) {
        Layer layer;
        if (layers.empty()) {
            layer.num_cols = input_size;
        } else {
            const Layer &last = layers.back();
            layer.num_cols = last.num_rows;
        }
        layers.emplace_back(std::move(layer));
        Layer &current = layers.back();

        current.num_rows = num_hidden;
        current.weights = std::make_unique<float[]>(current.num_rows * current.num_cols);
        current.biases = std::make_unique<float[]>(current.num_rows);
        current.pre_activ = std::make_unique<float[]>(current.num_rows);
        current.post_activ = std::make_unique<float[]>(current.num_rows);
        current.deltas = std::make_unique<float[]>(current.num_rows);
        //init everything but weights to zero
        std::fill(current.deltas.get(), current.deltas.get() + current.num_rows, 0.0f);
        std::fill(current.biases.get(), current.biases.get() + current.num_rows, 0.0f);

        const int rows = current.num_rows;

        float *in = current.pre_activ.get();
        float *out = current.post_activ.get();
        current.activation = [rows, in, out]() {
            Activ::apply(in, in + rows, out);
        };

        current.derivative = [rows](float *in, float *out) {
            Activ::apply_deriv(in, in + rows, out);
        };
    }

    void init_weights();

    template<typename Loss>
    void train(Train::Data &data, float lr, size_t epoch) {
        std::unique_ptr<float[]> input = std::make_unique<float[]>(input_size);
        std::unique_ptr<float[]> output = std::make_unique<float[]>(get_output_layer().num_rows);
        for (auto i = 0; i < epoch; ++i) {
            std::shuffle(data.mutable_data()->begin(), data.mutable_data()->end(), generator);

            for (auto k = 0; k < data.data_size(); ++k) {
                const Train::Point &point = data.data(k);
                for (auto p = 0; p < input_size; ++p) {
                    input[p] = point.input(p);
                }
                for (auto p = 0; p < get_output_layer().num_rows; ++p) {
                    output[p] = point.output(p);
                }
                forward_pass(input.get());

                back_propagate<Loss>(output.get());

                for (Layer &layer : layers) {
                    for (auto x = 0; x < layer.num_rows; ++x) {
                        layer.biases[x] -= lr * layer.deltas[x];
                    }
                }
                for (auto x = 0; x < layers.front().num_rows; ++x) {
                    for (auto y = 0; y < input_size; ++y) {

                        layers.front().weights[x * layers.front().num_cols + y] -=
                                lr * layers.front().deltas[x] * input[y];
                    }
                }

                for (auto l = 1; l < layers.size(); ++l) {
                    Layer &layer = layers[l];
                    Layer &previous = layers[l - 1];

                    for (auto x = 0; x < layer.num_rows; ++x) {
                        for (auto y = 0; y < layer.num_cols; ++y) {
                            layer.weights[x * layer.num_cols + y] -= lr * layer.deltas[x] * previous.post_activ[y];
                        }
                    }
                }

            }
        }
    }

    void forward_pass(float *input);

    Layer &get_output_layer();

    const Layer &get_output_layer() const;

    void train(Train::Data &data, float lr, size_t epoch);
};

#endif //NETWORK_NETWORK_H

max_stars_repo_name: CheckersGuy/Network

max_issues_repo_name: CheckersGuy/Network

max_forks_repo_name: CheckersGuy/Network


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_MATRIX_NORM_H
#define SCILIB_LINALG_MATRIX_NORM_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <lapacke.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <type_traits>

namespace Sci {
namespace Linalg {

namespace stdex = std::experimental;

// Matrix norm of a general rectangular matrix:
//
// Types of matrix norms:
// - M, m:       largest absolute value of the matrix
// - 1, O, o:    1-norm of the matrix (maximum column sum)
// - I, i:       infinity norm of the matrix (maximum row sum)
// - F, f, E, e: Frobenius norm of the matrix (square root of sum of squares)
//
template <class T, class Layout>
    requires(std::is_same_v<std::remove_cv_t<T>, double>)
inline auto matrix_norm(Sci::Matrix_view<T, Layout> a, char norm)
{
    static_assert(a.is_contiguous());

    assert(norm == 'M' || norm == 'm' || norm == '1' || norm == 'O' ||
           norm == 'o' || norm == 'I' || norm == 'i' || norm == 'F' ||
           norm == 'f' || norm == 'E' || norm == 'e');

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        lda = m;
    }
    return LAPACKE_dlange(matrix_layout, norm, m, n, a.data(), lda);
}

template <class Layout, class Allocator>
inline double matrix_norm(const Sci::Matrix<double, Layout, Allocator>& a,
                          char norm)
{
    return matrix_norm(a.view(), norm);
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_MATRIX_NORM_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
//
// Created by daeyun on 7/3/17.
//

#pragma once

#include <random>
#include <algorithm>
#include <gsl/gsl_assert>
#include <cstring>

namespace scene3d {
namespace Random {

static std::mt19937 &Engine() {
  thread_local static std::mt19937 engine{std::random_device{}()};
  return engine;
}

template<class Iter>
static void Shuffle(Iter begin, Iter end) {
  std::shuffle(begin, end, Engine());
}

static int UniformInt(int n) {
  std::uniform_int_distribution<decltype(n)> dist{0, n - 1};
  return dist(Engine());
}

static double Rand() {
  std::uniform_real_distribution<double> dist{0, 1};
  return dist(Engine());
}

// If n is greater than pool size, some items will be chosen more than once.
template<class T>
static void ChooseN(size_t n, std::vector<T> *mutable_pool, std::vector<T> *choices) {
  Expects(!mutable_pool->empty());
  auto left_size = mutable_pool->size();
  for (int i = 0; i < n; ++i) {
    std::string RandomString(size_t length);
    if (left_size <= 0) {
      left_size = mutable_pool->size();
    }
    auto &chosen = mutable_pool->at(UniformInt(left_size));
    choices->push_back(chosen);
    std::swap(chosen, mutable_pool->at(left_size - 1));
    --left_size;
  }
}

static std::string RandomString(size_t length) {
  auto randchar = []() -> char {
    constexpr char charset[] =
        "0123456789"
            "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
            "abcdefghijklmnopqrstuvwxyz";
    auto i = UniformInt(static_cast<int>(std::strlen(charset)));
    return charset[i];
  };
  std::string str(length, 0);
  std::generate_n(str.begin(), length, randchar);
  return str;
}

}
}

max_stars_repo_name: daeyun/mesh-point-cloud

max_issues_repo_name: daeyun/mesh-point-cloud

max_forks_repo_name: daeyun/mesh-point-cloud


text:
/*
 * File:   array_funcs.h
 * Author: fedro
 *
 * Created on December 10, 2014, 3:25 PM
 */

#ifndef ARRAY_FUNCS_H
#define	ARRAY_FUNCS_H

#include "array.h"
#include <assert.h>
#include <gsl/gsl_statistics_double.h>
#include <cmath>
#include <complex>
#include <fstream>

namespace neurophys {
namespace arr {

template<class T>
Array<T> zeros(const size_t size)
{
    // this assumes that whatever the type T, setting the memory byte-wise to
    // zero means the elements become 0

    Array<T> result(size);
    memset(result.data(), 0, size*sizeof(T));
    return result;
}

template<class T>
Array<T> ones(const size_t size)
{
    Array<T> result(size);
    for (int i = 0; i < size; i++)
        result[i] = 1.;
    return result;
}

template<class T>
const T sum(const Array<T>& ar)
{
    T res = 0.;
    for (int i = 0; i < ar.size(); i++)
        res += ar[i];
    return res;
}

template<class T>
const T max(const Array<T>& ar)
{
    T max = ar[0];
    for (int i = 0; i < ar.size(); i++)
        if (ar[i] > max)
            max = ar[i];
    return max;
}

template <class T>
void columns_to_stream(std::ostream& out, const Array<T>& x, const Array<T>& y)
{
    /*if (x.size() != y.size())
        throw ArraySizeMismatchException();*/
    for (int i = 0; i < std::min(x.size(),y.size()); i++)
        out << x[i] << "\t" << y[i] << "\n";
}

template <class T>
void columns_to_stream(std::ostream& out, const Array<T>& x, const Array<T>& y,
                        const Array<T>& z)
{
    /*if (x.size() != y.size())
        throw ArraySizeMismatchException();*/
    for (int i = 0; i < std::min(x.size(),y.size()); i++)
        out << x[i] << "\t" << y[i] << "\t" << z[i] << "\n";
}

template <class T>
void columns_to_file(const char filename[], const Array<T>& x, const Array<T>& y)
{
    std::ofstream of;
    of.open(filename, std::ios::out);
    columns_to_stream(of, x, y);
    of.close();
}

template <class T>
void columns_to_file(const char filename[], const Array<T>& x, const Array<T>& y,
                        const Array<T>& z)
{
    std::ofstream of;
    of.open(filename, std::ios::out);
    columns_to_stream(of, x, y, z);
    of.close();
}

template <class T>
Array<T> downsample(const Array<T>& src, const size_t destsize)
{
    assert(destsize <= src.size());
    Array<T> result = zeros<double>(destsize);
    const double f = static_cast<double>(destsize)/src.size();
    for (int i = 0; i < src.size(); i++)
    {
        result[i * f] += src[i] * f;
    }
    return result;
}

double mean(const Array<double>& ar);
double var(const Array<double>& ar);
Array<double> range(const double a, const double b, const double inc=1);
Array<double> cos(const Array<double>& x);
Array<double> sin(const Array<double>& x);
Array<std::complex<double> > conj(const Array<std::complex<double> >& x);
Array<double> real(const Array<std::complex<double> >& x);
Array<double> imag(const Array<std::complex<double> >& x);
Array<double> abs(const Array<std::complex<double> >& x);

}}

#endif	/* ARRAY_FUNCS_H */


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <cem.h>
#include <mpi.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

void randperm(int*, int);
void load_double_mat(double*, int, int, int, char*);
void load_int_mat(int*, int, int, int, char*);
void load_double_array(double*, int, char*);
void load_int_array(int*, int, char*);
void load_double_value(int,char);
void load_int_value(int,char);
int mat_size_row(int*);
int find_value_in_array(int*,int,int,int);
void mat2d_prod(double*, int, int, double*, int, int, double*, int, int);
double obtainerr2_par(int,int,double*, int, int, int*, int, double*);
void least_square_solver(double*, int, int, double*, double*);

#define MASTER 0
#define FROM_MASTER 1
#define FROM_WORKER 2

int main(int argc, char **argv)
{
    int n1, n2, n3;
    int np, nfu;
    double kv[3][3], pk[3][3], pp[8][3];
    double *rp, *rpL, *rpN, *rpO;
    int *map_to_cluster1, *map_to_cluster2, *map_to_cluster3;
    int *nlist;
    int neighbor_num;
    int ncluster, ncluster2, ncluster3, ncorr_col;
    int c2start, c3start;
    int data;
    double *E, *Ef;
    int ndata;
    int i, j, k, nj, nk, ctn;
    int iter, max_iter;
    double kT=0.0517;
    int howmanycluster;
    int dispfreq;

    int construct_corr_mat = 0;

    // mpi parameters and initialization
    int numprocs, rank, mtype;
    int row_dist_size;
    int row_ini, row_end, row_offset;
    MPI_Status status;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);  // Get # processors
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);      // Get my rank (id)

    if (numprocs>1)
        printf("Parallel computing: the input files are read by all nodes. Multiple messages are not due to error!\n",numprocs);

    // Construct conrrelation matrix. It can also be loaded.

    if (construct_corr_mat) {
        printf("Not activated yet.\n");
/*        // Task: make a function to load values from clusterlist run, later
        np = 48;            // for 3x2x2 cell
        nfu = np/2;         // for 3x2x2 cell
        neighbor_num = 47;  // for 3x2x2 cell
        ncluster2 = 10;
        ncluster3 = 62;

        // load lattice and cluster data
        load_double_mat(kv,3,3,1,"kv.dat");
        load_double_mat(rp,np,3,1,"rp.dat");
        load_double_mat(rpL,npL,3,1,"rpL.dat");
        load_double_mat(rpN,npN,3,1,"rpN.dat");
        load_double_mat(rpO,npO,3,1,"rpO.dat");
        load_double_array(cluster2,"cluster2.dat");
        load_double_array(cluster3,"cluster3.dat");
        load_int_array(map_to_cluster1,"map_to_cluster1.dat");
        load_int_array(map_to_cluster2,"map_to_cluster2.dat");
        load_int_array(map_to_cluster3,"map_to_cluster3.dat");
        load_int_array(nlist,"nlist.dat");

        ncluster  = ncluster2+ncluster3+1+1;
        ncorr_col = 1 + 2*np + 3*ncluster2 + 4*ncluster3;
        c2start = 3;    //2*length(cluster1)+1. Don't need -1: cluster2 idx starts from 0
        c3start = 3*ncluster2 + c2start;

        // load LNO occupation variable data
        load_int_array(data,"LNC_occu_data.dat");
        load_double_array(E,"E.dat");
        load_double_array(Ef,"Ef.dat");
        ndata = mat_size_row(data);
        int nL[ndata];
        int nC[ndata];
        int nN[ndata];

        for (i=0;i<ndata;i++) {
            nL[i] = find_value_in_array(data,i,0);
            nC[i] = find_value_in_array(data,i,1);
            nN[i] = find_value_in_array(data,i,-1);
        }

        // initialize the correlation matrix
        double corr_mat[ndata][ncorr_col];
        obtain_corr_mat(corr_mat, ndata, ncorr_col, data, c2start, c3start);

        // post data processing: excluding degenerate cases
        int nondegenerate[ndata];
        for (i=0;i<ndata;i++)
            nondegenerate[i] = 1;
        for (i=0;i<ndata;i++) {
            for (j=0;j<i;j++) {
                ctn = 1;
                for (l=0;l<ncorr_col;l++)
                    if (corr_mat[j][l]~=corr_mat[i][l]) {
                        ctn = 0;
                        break;
                    }
                if (ctn == 1)
                    if (Ef[j] < Ef[i])
                        nondegenerate[i] = 0;
                    else
                        nondegenerate[j] = 0;
            }
            for (j=i+1;j<ndata;j++) {
                ctn = 1;
                for (l=0;l<ncorr_col;l++)
                    if (corr_mat[j][l]~=corr_mat[i][l]) {
                        ctn = 0;
                        break;
                    }
                if (ctn == 1)
                    if (Ef[j] < Ef[i])
                        nondegenerate[i] = 0;
                    else
                        nondegenerate[j] = 0;
            }
        }

        int n_nondegenerate = 0;
        for (i=0;i<ndata;i++)
            if (nondegenerate[i] == 1)
                n_nondegenerate++;
        double corr_mat_u[n_nondegenerate][ncorr_col];
        double Ef_u[n_nondegenerate];
        double nC_r[n_nondegenerate];
        l = 0;
        for (i=0;i<ndata;i++)
            if (nondegenerate[i] == 1) {
                for (j=0;j<ncorr_col;j++)
                    corr_mat_u[l][j] = corr_mat[i][j];
                Ef_u[l] = Ef[i];
                nC_u[l] = nC[i];
                l++;
            }

        // select data sets near convex hull
        int n_corr_mat_ug_row = 0;
        int near_convh[n_nondegenerate];
        for (i=0;i<n_nondegenerate;i++)
            if (Ef_u[i] < ef_conv(nC[i]) + near_convh_cutoff) {
                near_convh[i] = 1;
                n_corr_mat_ug_row++;
            }
            else
                near_convh[i] = 0;

        double corr_mat_ug[n_corr_mat_ug_row][ncorr_col];
        double Ef_ug[n_corr_mat_ug_row];
        double nC_ug[n_corr_mat_ug_row];
        l = 0;
        for (i=0;i<n_corr_mat_ug_row)
            if (near_convh[i] == 1) {
                for (j=0;j<ncorr_col;j++)
                    corr_mat_ug[l][j] = corr_mat_u[i][j];
                Ef_ug[l] = Ef_u[i];
                nC_ug[l] = nC_u[i];
                l++;
            }

        // select the non-singular columns only
        int n_non_singular_col = 0;
        int col_singularity[ncorr_col];
        for (i=0;i<ncorr_col;i++) {
            ctn = 0;
            for (j=1;j<n_corr_mat_ug_row;j++)
                if (corr_mat_ug[j][i] ~= corr_mat_ug[0][i]) {
                    ctn = 1;
                    break;
                }
            if (ctn == 1) {
                col_singularity[i] = 1;  // 0:singular, 1:non-singular
                n_non_singular_col++;
            }
            else
                col_singularity[i] = 0;
        }
        double corr_mat_ugs[n_corr_mat_ug_row][n_non_singular_col];
        int non_singular_col_id[n_non_singular_col];
        l = 0;
        for (i=0;i<n_non_singular_col)
            if (col_singularity[i] == 1) {
                for (j=0;j<n_corr_mat_ug_row;j++)
                    corr_mat_ugs[j][l] = corr_mat_ug[j][l];
                non_singular_col_id[l] = i;
                l++;
            }
 */
    }
//    else {
    int n_corr_mat_ug_row;
    int n_non_singular_col;
    FILE *fp;
    char paramfilename[100]="param_findcluster.dat";
    char corr_mat_ugs_filename[200];
    fp = fopen(paramfilename, "r");
    fscanf(fp, "%s %d %d %d %d %lf %d",corr_mat_ugs_filename, &n_corr_mat_ug_row, &n_non_singular_col, &howmanycluster, &max_iter, &kT, &dispfreq);
    fclose(fp);
    printf("n_corr_mat_ug_row=%d, n_non_singular_col=%d, howmanycluster=%d, max_iter=%d, dispfreq=%d\n",n_corr_mat_ug_row, n_non_singular_col, howmanycluster, max_iter, dispfreq);
    double *corr_mat_ugs;
    corr_mat_ugs = (double*) malloc(n_corr_mat_ug_row*n_non_singular_col*sizeof(double));
    load_double_mat(corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,1,corr_mat_ugs_filename);

    int usefulcorr_col[n_non_singular_col];
    load_int_array(usefulcorr_col,n_non_singular_col,"usefulcorr_col.dat");

    double Ef_ug[n_corr_mat_ug_row];
    load_double_array(Ef_ug,n_corr_mat_ug_row,"Ef_ug.dat");

    int nC_ug[n_corr_mat_ug_row];
    load_int_array(nC_ug,n_corr_mat_ug_row,"nC_ug.dat");

//    }
    printf("step1\n");

    row_offset = n_corr_mat_ug_row/numprocs;
    double err2_sum_from_master;
    double err2_sum_from_worker;
    int cluster_set1[howmanycluster];

    /**** Start of MC simulation ****/
    if (rank == MASTER) {
        int cluster_set1_old[howmanycluster], cluster_set1_min[howmanycluster], cluster_set2[n_non_singular_col-howmanycluster];
        double cvs, cvs_old, cvs_min, cvs_tol;
        int select1, select2, target, candidate;

        cvs_tol = 0.01;
        for (i=0;i<howmanycluster;i++)
            cluster_set1[i] = i;
        for (i=0;i<n_non_singular_col-howmanycluster;i++)
            cluster_set2[i] = i+howmanycluster;

        mtype = FROM_MASTER;
        for (i=1;i<numprocs;i++) {
            MPI_Send(&cluster_set1, howmanycluster, MPI_INT, i, mtype, MPI_COMM_WORLD);
        }
        row_ini = rank * row_offset;
        row_end = row_ini +row_offset;
        err2_sum_from_master = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
        printf("initial err2 = %f\n",err2_sum_from_master);
        for (i=1;i<numprocs;i++) {
            mtype = FROM_WORKER;
            MPI_Recv(&err2_sum_from_worker, 1, MPI_DOUBLE, i, mtype, MPI_COMM_WORLD, &status);
            err2_sum_from_master += err2_sum_from_worker;
        }
        printf("initial err2 = %f\n",err2_sum_from_master);
        cvs = sqrt(err2_sum_from_master/n_corr_mat_ug_row);
        cvs_min = cvs;
        for (i=0;i<howmanycluster;i++)
            cluster_set1_min[i] = cluster_set1[i];
        printf("initial cvs = %f\n",cvs);
        for (iter=0;iter<max_iter;iter++) {
            select1 = (int) floor(drand48()*howmanycluster);
            select2 = (int) floor(drand48()*(n_non_singular_col-howmanycluster));
            target = cluster_set1[select1];
            candidate = cluster_set2[select2];
            cvs_old = cvs;
            cluster_set1[select1] = candidate;
            mtype = FROM_MASTER;
            for (i=1;i<numprocs;i++) {
                MPI_Send(&cluster_set1, howmanycluster, MPI_INT, i, mtype, MPI_COMM_WORLD);
            }
            row_ini = rank * row_offset;
            row_end = row_ini +row_offset;
            err2_sum_from_master = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
            for (i=1;i<numprocs;i++) {
                mtype = FROM_WORKER;
                MPI_Recv(&err2_sum_from_worker, 1, MPI_DOUBLE, i, mtype, MPI_COMM_WORLD, &status);
                err2_sum_from_master += err2_sum_from_worker;
            }
            cvs = sqrt(err2_sum_from_master/n_corr_mat_ug_row);
            if (cvs < cvs_old)
                cluster_set2[select2] = target;
            else if (exp(-(cvs-cvs_old)/kT) > drand48())
                cluster_set2[select2] = target;
            else {
                cluster_set1[select1] = target;
                cvs = cvs_old;
            }
            if (cvs < cvs_min) {
                for (i=0;i<howmanycluster;i++)
                    cluster_set1_min[i] = cluster_set1[i];
                cvs_min = cvs;
            }
            if (cvs < cvs_tol)
                break;
            if (iter%dispfreq==0)
                printf("iter=%8d, cvs=%f, cvs_min=%f, cvs_min/fu=%f\n",iter,cvs,cvs_min,cvs_min/72);
        }

        // print the result
        printf("End of MC iteration. The selected clusters are...\n");
        printf("ID in non-singular | ID in entire\n");
        for (i=0;i<howmanycluster;i++)
            printf("    %4d                %4d\n",cluster_set1_min[i],usefulcorr_col[cluster_set1_min[i]]);

        double err_pred, rms;
        double Ef_pred[n_corr_mat_ug_row];
        double *corr_mat_ugs_reduced;
        corr_mat_ugs_reduced = (double*) malloc(n_corr_mat_ug_row*howmanycluster*sizeof(double));
        double eci[howmanycluster];
        for (i=0;i<n_corr_mat_ug_row;i++) {
            for (j=0;j<howmanycluster;j++)
                *(corr_mat_ugs_reduced+howmanycluster*i+j) = *(corr_mat_ugs+n_non_singular_col*i+cluster_set1_min[j]);
        }
        least_square_solver(corr_mat_ugs_reduced,n_corr_mat_ug_row,howmanycluster,Ef_ug,eci);
        mat2d_prod(Ef_pred,n_corr_mat_ug_row,1,corr_mat_ugs_reduced,n_corr_mat_ug_row,howmanycluster,eci,howmanycluster,1);
        err_pred = 0.0;
        for (i=0;i<n_corr_mat_ug_row;i++)
            err_pred += (Ef_ug[i] - Ef_pred[i])*(Ef_ug[i] - Ef_pred[i]);
        rms = sqrt(err_pred/n_corr_mat_ug_row);

        // Write result into files
        FILE *fp2;
        char cluster_set_filename[200];
        sprintf(cluster_set_filename, "cluster_set_magN%d.dat",howmanycluster);
        fp2 = fopen(cluster_set_filename,"w");
        for (i=0;i<howmanycluster;i++)
            fprintf(fp2,"%4d\t%4d\t%.4f\n",cluster_set1_min[i],usefulcorr_col[cluster_set1_min[i]],eci[i]);
    }
    else {
        for (iter=0;iter<(max_iter+1);iter++) {
            mtype = FROM_MASTER;
            MPI_Recv(&cluster_set1, howmanycluster, MPI_INT, MASTER, mtype, MPI_COMM_WORLD, &status);
            row_ini = rank * row_offset;
            row_end = row_ini + row_offset;
            err2_sum_from_worker = obtainerr2_par(row_ini,row_end,corr_mat_ugs,n_corr_mat_ug_row,n_non_singular_col,cluster_set1,howmanycluster,Ef_ug);
            mtype = FROM_WORKER;
            MPI_Send(&err2_sum_from_worker, 1, MPI_DOUBLE, MASTER, mtype, MPI_COMM_WORLD);
        }
    }
    MPI_Finalize();

}

void load_double_mat(double *A, int Arow, int Acol, int Apgs, char *datfilename) {

    FILE *fp;
    int i, j, k;
    double tmp;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    for (k=0;k<Apgs;k++)
        for (i=0;i<Arow;i++)
            for (j=0;j<Acol;j++) {
                fscanf(fp, "%lf", &tmp);
                *(A+Arow*Acol*k+Acol*i+j) = tmp;
            }

    fclose(fp);

}

void load_int_mat(int *A, int Arow, int Acol, int Apgs, char *datfilename) {

    FILE *fp;
    int i, j, k, tmp;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    for (k=0;k<Apgs;k++)
        for (i=0;i<Arow;i++)
            for (j=0;j<Acol;j++) {
                fscanf(fp, "%d", &tmp);
                *(A+Arow*Acol*k+Acol*i+j) = tmp;
            }

    fclose(fp);

}

void load_int_array(int *array, int arraysize, char *datfilename) {

    FILE *fp;
    int i;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    //    printf("step2\n");
    for (i=0;i<arraysize;i++) {
        fscanf(fp, "%d", &array[i]);
        //        printf("%d\n",array[i]);
    }

    fclose(fp);
    //    printf("step3\n");
}

void load_double_array(double *array, int arraysize, char *datfilename) {

    FILE *fp;
    int i;

    printf("filename: %s\n", datfilename);
    fp = fopen(datfilename, "r");
    //    printf("step4\n");
    for (i=0;i<arraysize;i++) {
        fscanf(fp, "%lf", &array[i]);
        //        printf("%f\n",array[i]);
    }

    fclose(fp);
    //    printf("step5\n");
}

void randperm(int *a, int np) {

    int i, id1, id2, id3, tmpr, num_perm;

    num_perm = 100*np;

    for (i=0;i<np;i++)
        *(a+i) = i;

    for (i=0;i<num_perm;i++) {
        id1 = rand() % np;
        id3 = rand() % np;
        id2 = (id1 + id3) % np;
        if (id1 != id2) {
            tmpr = *(a+id2);
            *(a+id2) = *(a+id1);
            *(a+id1) = tmpr;
        }
    }
}

int find_value_in_array(int *A, int target_rowA, int colA, int target_value) {
    int i;

    for (i=0;i<colA;i++)
        if (*(A+colA*target_rowA+i)==target_value)
            return(i);
    if (i==colA)
        return(-1);
}


text:
#pragma once
#include "../utils/utils.h"
#include <gsl/gsl>
#include <limits>

namespace Halley {
	class Compression {
	public:
		static Bytes compress(const Bytes& bytes);
		static Bytes compress(gsl::span<const gsl::byte> bytes);
		static Bytes decompress(const Bytes& bytes, size_t maxSize = std::numeric_limits<size_t>::max());
		static Bytes decompress(gsl::span<const gsl::byte> bytes, size_t maxSize = std::numeric_limits<size_t>::max());
		static std::shared_ptr<const char> decompressToSharedPtr(gsl::span<const gsl::byte> bytes, size_t& outSize, size_t maxSize = std::numeric_limits<size_t>::max());

		static Bytes compressRaw(gsl::span<const gsl::byte> bytes, bool insertLength);
		static Bytes decompressRaw(gsl::span<const gsl::byte> bytes, size_t maxSize, size_t expectedSize = 0);
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley