text:
// Copyright (c) 2013-2016 Anton Kozhevnikov, Thomas Schulthess
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that
// the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
//    following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
//    and the following disclaimer in the documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

/** \file sirius_internal.h
 *
 *  \brief Contains basic definitions and declarations.
 */

#ifndef __SIRIUS_INTERNAL_H__
#define __SIRIUS_INTERNAL_H__

#include <omp.h>
#include <signal.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf_bessel.h>
#include <fftw3.h>
#include <vector>
#include <complex>
#include <iostream>
#include <algorithm>
#include "config.h"
#include "communicator.hpp"
#include "gpu.h"
#include "runtime.h"

#ifdef __PLASMA
extern "C" void plasma_init(int num_cores);
#endif

#ifdef __LIBSCI_ACC
extern "C" void libsci_acc_init();
extern "C" void libsci_acc_finalize();
#endif

/// Namespace of the SIRIUS library.
namespace sirius {

    inline void initialize(bool call_mpi_init__)
    {
        if (call_mpi_init__) Communicator::initialize();

        #ifdef __GPU
        cuda_create_streams(omp_get_max_threads() + 1);
        cublas_create_handles(omp_get_max_threads() + 1);
        #endif
        #ifdef __MAGMA
        magma_init_wrapper();
        #endif
        #ifdef __PLASMA
        plasma_init(omp_get_max_threads());
        #endif
        #ifdef __LIBSCI_ACC
        libsci_acc_init();
        #endif

        assert(sizeof(int) == 4);
        assert(sizeof(double) == 8);
    }

    inline void finalize()
    {
        Communicator::finalize();
        #ifdef __MAGMA
        magma_finalize_wrapper();
        #endif
        #ifdef __LIBSCI_ACC
        libsci_acc_finalize();
        #endif
        #ifdef __GPU
        cublas_destroy_handles(omp_get_max_threads() + 1);
        cuda_destroy_streams();
        cuda_device_reset();
        #endif
        fftw_cleanup();
    }
};

#define TERMINATE_NO_GPU TERMINATE("not compiled with GPU support");

#define TERMINATE_NO_SCALAPACK TERMINATE("not compiled with ScaLAPACK support");

#define TERMINATE_NOT_IMPLEMENTED TERMINATE("feature is not implemented");

#endif // __SIRIUS_INTERNAL_H__

/** \mainpage Welcome to SIRIUS
 *  \section intro Introduction
 *  SIRIUS is a domain-specific library for electronic structure calculations. It supports full-potential linearized
 *  augmented plane wave (FP-LAPW) and pseudopotential plane wave (PP-PW) methods and is designed to work with codes
 *  such as Exciting, Elk, Quantum ESPRESSO, etc.
 *  \section install Installation
 *   ...
 */

/** \page stdvarname Standard variable names
 *
 *  Below is the list of standard names for some of the loop variables:
 *
 *  l - index of orbital quantum number \n
 *  m - index of azimutal quantum nuber \n
 *  lm - combined index of (l,m) quantum numbers \n
 *  ia - index of atom \n
 *  ic - index of atom class \n
 *  iat - index of atom type \n
 *  ir - index of r-point \n
 *  ig - index of G-vector \n
 *  idxlo - index of local orbital \n
 *  idxrf - index of radial function \n
 *  xi - compbined index of lm and idxrf (product of angular and radial functions) \n
 *  ik - index of k-point \n
 *  itp - index of (theta, phi) spherical angles \n
 *
 *  The _loc suffix is often added to the variables to indicate that they represent the local fraction of the elements
 *  assigned to the given MPI rank.
 */

//! \page coding Coding style
//!
//! Below are some basic style rules that we follow:
//!     - Page width is approximately 120 characters. Screens are wide nowdays and 80 characters is an
//!       obsolete restriction. Going slightly over 120 characters is allowed if it is requird for the line continuity.
//!     - Identation: 4 spaces (no tabs)
//!     - Coments are inserted before the code with slash-star style starting with the lower case:
//!       \code{.cpp}
//!           /* call a very important function */
//!           do_something();
//!       \endcode
//!     - Spaces between most operators:
//!       \code{.cpp}
//!           if (i < 5) {
//!               j = 5;
//!           }
//!
//!           for (int k = 0; k < 3; k++)
//!
//!           int lm = l * l + l + m;
//!
//!           double d = std::abs(e);
//!
//!           int k = idx[3];
//!       \endcode
//!     - Spaces between function arguments:
//!       \code{.cpp}
//!           double d = some_func(a, b, c);
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           double d=some_func(a,b,c);
//!       \endcode
//!       or
//!       \code{.cpp}
//!           double d = some_func( a, b, c );
//!       \endcode
//!     - Spaces between template arguments, but not between <> brackets:
//!       \code{.cpp}
//!           std::vector<std::array<int, 2>> vec;
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           std::vector< std::array< int, 2 > > vec;
//!       \endcode
//!     - Curly braces for classes and functions start form the new line:
//!       \code{.cpp}
//!           class A
//!           {
//!               ....
//!           };
//!
//!           inline int num_points()
//!           {
//!               return num_points_;
//!           }
//!       \endcode
//!     - Curly braces for if-statements, for-loops, switch-case statements, etc. start at the end of the line:
//!       \code{.cpp}
//!           for (int i: {0, 1, 2}) {
//!               some_func(i);
//!           }
//!
//!           if (a == 0) {
//!               printf("a is zero");
//!           } else {
//!               printf("a is not zero");
//!           }
//!
//!           switch (i) {
//!               case 1: {
//!                   do_something();
//!                   break;
//!               case 2: {
//!                   do_something_else();
//!                   break;
//!               }
//!           }
//!       \endcode
//!     - Even single line 'if' statements and 'for' loops must have the curly brackes:
//!       \code{.cpp}
//!           if (i == 4) {
//!               some_variable = 5;
//!           }
//!
//!           for (int k = 0; k < 10; k++) {
//!               do_something(k);
//!           }
//!       \endcode
//!     - Reference and pointer symbols are part of type:
//!       \code{.cpp}
//!           std::vector<double>& vec = make_vector();
//!
//!           double* ptr = &vec[0];
//!
//!           auto& atom = unit_cell().atom(ia);
//!       \endcode
//!     - Const modifier follows the type declaration:
//!       \code{.cpp}
//!           std::vector<int> const& idx() const
//!           {
//!               return idx_;
//!           }
//!       \endcode
//!     - Names of class members end with underscore:
//!       \code{.cpp}
//!           class A
//!           {
//!               private:
//!                   int lmax_;
//!           };
//!       \endcode
//!     - Setter method starts from set_, getter method is a variable name itself:
//!       \code{.cpp}
//!           class A
//!           {
//!               private:
//!                   int lmax_;
//!               public:
//!                   int lmax() const
//!                   {
//!                       return lmax_;
//!                   }
//!                   void set_lmax(int lmax__)
//!                   {
//!                       lmax_ = lmax__;
//!                   }
//!           };
//!       \endcode
//!     - Single-line functions should not be flattened:
//!       \code{.cpp}
//!           struct A
//!           {
//!               int lmax() const
//!               {
//!                   return lmax_;
//!               }
//!           };
//!       \endcode
//!       but not
//!       \code{.cpp}
//!           struct A
//!           {
//!               int lmax() const { return lmax_; }
//!           };
//!       \endcode
//!     - Header guards have a standard name: double underscore + file name in capital letters + double underscore
//!       \code{.cpp}
//!           #ifndef __SIRIUS_INTERNAL_H__
//!           #define __SIRIUS_INTERNAL_H__
//!           ...
//!           #endif // __SIRIUS_INTERNAL_H__
//!       \endcode
//! We use clang-format utility to enforce the basic formatting style. Please have a look at .clang-format config file
//! in the source root folder for the definitions.
//!
//! Class naming convention.
//!
//! Problem: all 'standard' naming conventions are not satisfactory. For example, we have a class
//! which does a DFT ground state. Following the common naming conventions it could be named like this:
//! DFTGroundState, DftGroundState, dft_ground_state. Last two are bad, because DFT (and not Dft or dft)
//! is a well recognized abbreviation. First one is band because capital G adds to DFT and we automaticaly
//! read DFTG round state.
//!
//! Solution: we can propose the following: DFTgroundState or DFT_ground_state. The first variant still
//! doens't look very good because one of the words is captalized (State) and one (ground) - is not. So we pick
//! the second variant: DFT_ground_state (by the way, this is close to the Bjarne Stroustrup's naiming convention,
//! where he uses first capital letter and underscores, for example class Io_obj).
//!
//! Some other examples:
//!     - class Ground_state (composed of two words)
//!     - class FFT_interface (composed of an abbreviation and a word)
//!     - class Interface_XC (composed of a word and abbreviation)
//!     - class Spline (single word)
//!
//! Exceptions are allowed if it makes sense. For example, low level utility classes like 'mdarray' (multi-dimentional
//! array) or 'pstdout' (parallel standard output) are named with small letters.
//!

ext: h

max_stars_repo_name: cocteautwins/SIRIUS-develop

max_issues_repo_name: cocteautwins/SIRIUS-develop

max_forks_repo_name: cocteautwins/SIRIUS-develop


text:
#pragma once

#include "file_descriptor.h"
#include <filesystem>
#include <gsl/span>
#include <random>

namespace dogbox
{
    inline void create_random_file(std::filesystem::path const &file, uint64_t const size)
    {
        file_descriptor const created = create_file(file).value();
        file_descriptor const random = open_file_for_reading("/dev/urandom").value();
        std::array<std::byte, 0x10000> buffer;
        uint64_t written = 0;
        while (written < size)
        {
            size_t const reading = static_cast<size_t>(std::min<uint64_t>(buffer.size(), (size - written)));
            ssize_t const read_result = read(random.handle, buffer.data(), reading);
            if (read_result < 0)
            {
                TO_DO();
            }
            ssize_t const write_result = write(created.handle, buffer.data(), reading);
            if (write_result < 0)
            {
                TO_DO();
            }
            written += reading;
        }
    }
}

max_stars_repo_name: TyRoXx/dogbox

max_issues_repo_name: TyRoXx/dogbox

max_forks_repo_name: TyRoXx/dogbox


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>

#include "dnest.h"
#include "model3.h"

int which_level_update;
int num_params;
DNestFptrSet *fptrset_thismodel3;

void model3()
{
  int i, argc=0, narg=6;
  char **argv;

  argv = malloc(narg*sizeof(char *));
  for(i=0; i<narg; i++)
  {
    argv[i] = malloc(200*sizeof(char));
  }
  strcpy(argv[argc++], "dnest");
  strcpy(argv[argc++], "-s");
  strcpy(argv[argc++], "restart_dnest3.txt");
  strcpy(argv[argc++], "-l"); //level-dependnet sampling
  strcpy(argv[argc++], "-g"); // tag
  strcpy(argv[argc++], "3");

  /* setup szie of modeltype, which is used for dnest */
  num_params = 2;
  fptrset_thismodel3 = dnest_malloc_fptrset();

  /* setup functions used for dnest*/
  fptrset_thismodel3->from_prior = from_prior_thismodel3;
  fptrset_thismodel3->log_likelihoods_cal = log_likelihoods_cal_thismodel3;
  fptrset_thismodel3->log_likelihoods_cal_initial = log_likelihoods_cal_thismodel3;
  fptrset_thismodel3->log_likelihoods_cal_restart = log_likelihoods_cal_thismodel3;
  fptrset_thismodel3->perturb = perturb_thismodel3;
  fptrset_thismodel3->print_particle = print_particle_thismodel3;
  fptrset_thismodel3->restart_action = restart_action_model3;

  /* run dnest */
  dnest(argc, argv, fptrset_thismodel3, num_params, NULL, NULL, NULL, "./", "OPTIONS3", NULL, NULL);

  /* free memory */
  dnest_free_fptrset(fptrset_thismodel3);

  for(i=0; i<narg; i++)
    free(argv[i]);
  free(argv);

}

/*====================================================*/
/* users responsible for following struct definitions */

void from_prior_thismodel3(void *model)
{
  int i;
  double *params = (double *)model;
  for(i=0; i<num_params; i++)
  {
    params[i] = -6.0 + 12.0 * dnest_rand();
    dnest_wrap(&params[i], -6.0, 6.0);
  }
}

double log_likelihoods_cal_thismodel3(const void *model)
{
  double *params = (double *)model;
  double logL;
  double logl1;
  double logl2;

  logl1 = -0.5 * pow( (sqrt((params[0] - 3.0)* (params[0] - 3.0) + params[1]*params[1])- 2.0)/0.1,  2.0)
          - 0.5 * log(2.0*M_PI * 0.01);
  logl2 = -0.5 * pow( (sqrt((params[0] + 3.0)* (params[0] + 3.0) + params[1]*params[1])- 2.0)/0.1,  2.0)
          - 0.5 * log(2.0*M_PI * 0.01);

  double max = fmax(logl1, logl2);

  logL = log( exp(logl1-max) + exp(logl2-max) ) + max;

  return logL;
}

double perturb_thismodel3(void *model)
{
  double *params = (double *)model;
  double logH = 0.0, width, limit1, limit2;
  int which = dnest_rand_int(num_params), which_level;
  int size_levels;

  which_level_update = dnest_get_which_level_update();
  size_levels = dnest_get_size_levels();
  which_level = which_level_update > (size_levels - 20)?(size_levels-20):which_level_update;
  which_level = 0;
  if(which_level > 0 )
  {
    limit1 = limits[(which_level_update-1) * num_params *2 + which *2 ];
    limit2 = limits[(which_level_update-1) * num_params *2 + which *2 + 1];
  }
  else
  {
    limit1 = -6.0;
    limit2 = 6.0;
  }
  width = (limit2 - limit1);
  params[which] += width * dnest_randh();
  dnest_wrap(&params[which], -6.0, 6.0);
  return logH;
}
/*=======================================================*/

void print_particle_thismodel3(FILE *fp, const void *model)
{
  int i;
  double *params = (double *)model;
  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%f ", params[i]);
  }
  fprintf(fp, "\n");
}

/*========================================================*/

void restart_action_model3(int iflag)
{
  return;
}

ext: c

max_stars_repo_name: LiyrAstroph/CDNest


text:
#pragma once

#include <gsl.h>

#include <mitkCommon.h>
#include <mitkBaseData.h>

#include <vtkSmartPointer.h>
#include <vtkDataArray.h>

#include "SolverSetupServiceExports.h"

namespace crimson
{

/*! \brief   A class used for storing the simulation results and materials. */
class SolverSetupService_EXPORT SolutionData : public mitk::BaseData
{
public:
    mitkClassMacro(SolutionData, mitk::BaseData);
    mitkNewMacro1Param(Self, vtkSmartPointer<vtkDataArray>);
    //     itkCloneMacro(Self);
    //     mitkCloneMacro(Self);

    /*!
     * \brief   Gets the data in form of a vtkDataArray.
     */
    vtkDataArray* getArrayData() const;

protected:
    SolutionData(vtkSmartPointer<vtkDataArray> data);
    virtual ~SolutionData() {}

    SolutionData(const Self& other) = delete;

private:
    vtkSmartPointer<vtkDataArray> _data;

    void SetRequestedRegion(const itk::DataObject*) override {}
    void SetRequestedRegionToLargestPossibleRegion() override {}
    bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return true; }
    bool VerifyRequestedRegion() override { return true; }
};
}

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
#include <asf_terrcorr.h>
#include <stdio.h>
#include <assert.h>
#include <gsl/gsl_spline.h>
#include <asf_raster.h>

static float *
read_dem(meta_parameters *meta_dem, const char *demImg)
{
  int ns = meta_dem->general->sample_count;
  int nl = meta_dem->general->line_count;
  float *demData = MALLOC(sizeof(float)*ns*nl);

  FILE *fp = FOPEN(demImg, "rb");

  int ii;
  for (ii=0; ii<nl; ++ii) {
    get_float_line(fp, meta_dem, ii, demData + ii*ns);
    asfLineMeter(ii,nl);
  }

  FCLOSE(fp);
  return demData;
}

static void
sar_to_dem(meta_parameters *meta_sar, meta_parameters *meta_dem,
           double line_sar, double samp_sar,
           double *line_dem, double *samp_dem)
{
  double lat, lon;
  int ret;

  ret = meta_get_latLon(meta_sar, line_sar, samp_sar, 0, &lat, &lon);
  if (ret != 0) {
    asfPrintError("meta_get_latLon error: line = %f, samp = %f\n",
                  line_sar, samp_sar);
  }
  ret = meta_get_lineSamp(meta_dem, lat, lon, 0, line_dem, samp_dem);
  if (ret != 0) {
    asfPrintError("meta_get_lineSamp error: lat = %f, lon = %f\n",
                  lat, lon);
  }
}

static double
bilinear_interp_fn(double y, double x,
                   double p00, double p10, double p01, double p11)
{
  double x1 = 1-x;
  double y1 = 1-y;

  return p00*x1*y1 + p10*x*y1 + p01*x1*y + p11*x*y;
}

static float
interp_demData(float *demData, int nl, int ns, double l, double s)
{
  if (l<0 || l>=nl-1 || s<0 || s>=ns-1) {
    return 0;
  }

  int ix = (int)s;
  int iy = (int)l;

  int bilinear = l<3 || l>=nl-3 || s<3 || s>=ns-3;
  //int bilinear = 1;
  if (bilinear) {

    float p00 = demData[ix   + ns*(iy  )];
    float p10 = demData[ix+1 + ns*(iy  )];
    float p01 = demData[ix   + ns*(iy+1)];
    float p11 = demData[ix+1 + ns*(iy+1)];

    return (float)bilinear_interp_fn(s-ix, l-iy, p00, p10, p01, p11);
  }
  else {

    double ret, x[4], y[4], xi[4], yi[4];
    int ii;

    for (ii=0; ii<4; ++ii) {
      y[0] = demData[ix-1 + ns*(iy+ii-1)];
      y[1] = demData[ix   + ns*(iy+ii-1)];
      y[2] = demData[ix+1 + ns*(iy+ii-1)];
      y[3] = demData[ix+2 + ns*(iy+ii-1)];

      x[0] = ix - 1;
      x[1] = ix;
      x[2] = ix + 1;
      x[3] = ix + 2;

      gsl_interp_accel *acc = gsl_interp_accel_alloc ();
      gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, 4);
      gsl_spline_init (spline, x, y, 4);
      yi[ii] = gsl_spline_eval(spline, s, acc);
      gsl_spline_free (spline);
      gsl_interp_accel_free (acc);

      xi[ii] = iy + ii - 1;
    }

    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, 4);
    gsl_spline_init (spline, xi, yi, 4);
    ret = gsl_spline_eval(spline, l, acc);
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);

    return (float)ret;
  }

  asfPrintError("Impossible.");
}

static void check(int num, double a, double b)
{
  if (a==b || fabs(a-b)<0.0000001)
    return;

  asfPrintError("Failed test %d -- %f %f\n", num, a, b);
}

static void test_interp()
{
  check(1, bilinear_interp_fn(.5,.5,1,1,1,1), 1);
  check(2, bilinear_interp_fn(.1,.1,1,1,1,1), 1);
  check(3, bilinear_interp_fn(.5,.5,1,1,2,2), 1.5);
  check(4, bilinear_interp_fn(.5,.5,1,2,1,2), 1.5);
  check(5, bilinear_interp_fn(.5,.5,1,2,1,2), 1.5);
  check(6, bilinear_interp_fn(0,0,1,2,3,4), 1);
  check(7, bilinear_interp_fn(1,1,1,2,3,4), 4);
  check(8, bilinear_interp_fn(0,1,1,2,3,4), 2);
  check(9, bilinear_interp_fn(1,0,1,2,3,4), 3);
}

static void xy_interp(int line, int samp,
                      int line_lo, int line_hi, int samp_lo, int samp_hi,
                      double *lines, double *samps,
                      double *interp_line, double *interp_samp)
{
  assert(samp>=samp_lo && samp<=samp_hi);
  assert(line>=line_lo && line<=line_hi);

  double l = (double)(line - line_lo) / (double)(line_hi - line_lo);
  double s = (double)(samp - samp_lo) / (double)(samp_hi - samp_lo);

  *interp_line = bilinear_interp_fn(l, s, lines[0], lines[1], lines[2], lines[3]);
  *interp_samp = bilinear_interp_fn(l, s, samps[0], samps[1], samps[2], samps[3]);
}

static double
calc_err(meta_parameters *meta_sar, meta_parameters *meta_dem,
         int line_lo, int line_hi, int samp_lo, int samp_hi,
         double *lines, double *samps)
{
  int l = (line_lo + line_hi) * .5;
  int s = (samp_lo + samp_hi) * .5;

  double interp_line, interp_samp, real_line, real_samp;
  xy_interp(l, s, line_lo, line_hi, samp_lo, samp_hi, lines, samps,
            &interp_line, &interp_samp);

  sar_to_dem(meta_sar, meta_dem, l, s, &real_line, &real_samp);

  asfPrintStatus("  Interp: (%f,%f)\n", interp_line, interp_samp);
  asfPrintStatus("  Actual: (%f,%f)\n", real_line, real_samp);

  return hypot(interp_line - real_line, interp_samp - real_samp);
}

static void
get_interp_params(meta_parameters *meta_sar, meta_parameters *meta_dem,
                  int line_lo, int line_hi, int samp_lo, int samp_hi,
                  double *lines, double *samps)
{
  sar_to_dem(meta_sar, meta_dem, line_lo, samp_lo, &lines[0], &samps[0]);
  sar_to_dem(meta_sar, meta_dem, line_lo, samp_hi, &lines[1], &samps[1]);
  sar_to_dem(meta_sar, meta_dem, line_hi, samp_lo, &lines[2], &samps[2]);
  sar_to_dem(meta_sar, meta_dem, line_hi, samp_hi, &lines[3], &samps[3]);
}

static int
find_grid_size(meta_parameters *meta_sar, meta_parameters *meta_dem,
               int initial_size, double tolerance)
{
  // pick a square centered on the center of the image
  // we will start with the biggest square and compare the error at the center
  // point between using (1) bilinear interp and (2) the real mapping
  // if the error is larger than the tolerance, shrink square and try again
  int line = meta_sar->general->line_count / 2;
  int samp = meta_sar->general->sample_count / 2;
  int sz = initial_size;
  double lines[4], samps[4];

  while (1) {
    int sz2 = sz/2;
    int line_lo = line - sz2;
    int line_hi = line_lo + sz;
    int samp_lo = samp - sz2;
    int samp_hi = samp_lo + sz;
    get_interp_params(meta_sar, meta_dem, line_lo, line_hi, samp_lo, samp_hi,
                      lines, samps);
    double err = calc_err(meta_sar, meta_dem, line_lo, line_hi, samp_lo, samp_hi,
                          lines, samps);
    if (err < tolerance) {
      asfPrintStatus("Using square size %d (err %f < %f)\n", sz, err, tolerance);
      break;
    }
    asfPrintStatus("Square size %d, no good (err %f > %f)\n", sz, err, tolerance);
    sz = sz2;
    if (sz <= 2) {
      asfPrintStatus("Proceeding with grid size 1.\n");
      sz = 1;
      break;
    }
  }
  return sz;
}

int make_gr_dem(meta_parameters *meta_sar, const char *demBase, const char *output_name)
{
  char *demImg = appendExt(demBase, ".img");
  char *demMeta = appendExt(demBase, ".meta");
  int ret = make_gr_dem_ext(meta_sar, demImg, demMeta, 0, .1, output_name, 0);
  FREE(demImg);
  FREE(demMeta);
  return ret;
}

// This is the external facing function from this file.
// Given the geometry from a sar image (meta_sar), and a DEM (in the files
// demImg and demMeta), we write out an image/meta pair with the specified
// output_name, which is a subset of the given DEM that covers the given
// sar geometry.  This function uses interpolations for speed, so there is
// a tolerance parameter that specifies how accurate the line/sample mapping
// from sar to dem geometry needs to be.

// Input:
//   meta_parameters *meta_sar--  SAR geometry to subset the DEM
//   const char *demImg       --  DEM data filename
//   const char *demMeta      --  DEM metadata filename
//   int pad                  --  number of lines to add at the top/bottom/left/right
//   double tolerance         --  how accurate the approximation mapping needs to be,
//                                in units of pixels
//   const char *output_name  --  output filename (basename)
//   int test_mode            --  adds checks for the accuracy of the mapping, and
//                                does some unit testing
// Output:
//   no output parameters, the output is the output_name files (.img and .meta)
// Return Value:
//   return TRUE on success, FALSE on fail
//
int make_gr_dem_ext(meta_parameters *meta_sar, const char *demImg, const char *demMeta,
                    int pad, double tolerance, const char *output_name, int test_mode)
{
  if (test_mode)
    test_interp();

  asfPrintStatus("Reading DEM...\n");
  meta_parameters *meta_dem = meta_read(demMeta);
  float *demData = read_dem(meta_dem, demImg);
  int dnl = meta_dem->general->line_count;
  int dns = meta_dem->general->sample_count;

  char *outImg = appendExt(output_name, ".img");
  char *output_name_tmp, *outImgTmp;

  // do not do DEM smoothing if the DEM pixel size is better or close to the
  // SAR image's pixel size.
  int do_averaging = TRUE;
  if (meta_dem->general->y_pixel_size - 10 < meta_sar->general->y_pixel_size)
    do_averaging = FALSE;
  asfPrintStatus("Averaging: %s (DEM %f, SAR: %f)\n", do_averaging ? "YES" : "NO",
                 meta_dem->general->y_pixel_size,
                 meta_sar->general->y_pixel_size);
  if (do_averaging) {
    output_name_tmp = appendStr(output_name, "_unsmoothed");
    outImgTmp = appendExt(output_name_tmp, ".img");
  }
  else {
    output_name_tmp = STRDUP(output_name);
    outImgTmp = STRDUP(outImg);
  }

  // add the padding if requested
  meta_parameters *meta_out = meta_copy(meta_sar);
  meta_out->general->line_count += pad*2;
  meta_out->general->sample_count += pad*2;
  meta_out->general->start_line -= pad;
  meta_out->general->start_sample -= pad;

  // fixing up the output metadata.  Note that we must keep the SAR section
  // intact since that specifies our geometry which is the whole point of
  // this exercise.
  strcpy(meta_out->general->basename, meta_dem->general->basename);
  strcpy(meta_out->general->sensor, MAGIC_UNSET_STRING);
  strcpy(meta_out->general->processor, MAGIC_UNSET_STRING);
  strcpy(meta_out->general->mode, MAGIC_UNSET_STRING);
  strcpy(meta_out->general->sensor_name, MAGIC_UNSET_STRING);
  meta_out->general->image_data_type = DEM;
  meta_out->general->radiometry = MAGIC_UNSET_INT;
  strcpy(meta_out->general->acquisition_date, meta_dem->general->acquisition_date);
  meta_out->general->orbit = MAGIC_UNSET_INT;
  meta_out->general->orbit_direction = MAGIC_UNSET_CHAR;
  meta_out->general->frame = MAGIC_UNSET_INT;
  meta_out->general->band_count = 1;
  strcpy(meta_out->general->bands, "DEM");

  int nl = meta_out->general->line_count;
  int ns = meta_out->general->sample_count;

  // finding the right grid size
  int size = find_grid_size(meta_sar, meta_dem, 512, .1*tolerance);

  asfPrintStatus("Creating ground range image...\n");

  float *buf = MALLOC(sizeof(float)*ns*size);
  FILE *fpOut = FOPEN(outImgTmp, "wb");

  // these are for tracking the quality of the bilinear interp
  // not used if test_mode is false
  int num_out_of_tol = 0;
  int num_checked = 0;
  int num_bad = 0;
  double max_err = 0;
  double avg_err = 0;

  int ii, jj;
  for (ii=0; ii<nl; ii += size) {
    int line_lo = ii;
    int line_hi = ii + size;

    for (jj=0; jj<ns; jj += size) {
      double lines[4], samps[4];

      int samp_lo = jj;
      int samp_hi = jj + size;

      get_interp_params(meta_sar, meta_dem, line_lo, line_hi, samp_lo, samp_hi,
                        lines, samps);

      int iii, jjj;
      for (iii=0; iii<size; ++iii) {
        for (jjj=0; jjj<size && jj+jjj<ns; ++jjj) {
          int index = iii*ns + jj + jjj;
          assert(index < ns*size);

          double line_out, samp_out;
          xy_interp(ii+iii, jj+jjj, line_lo, line_hi, samp_lo, samp_hi, lines, samps,
                    &line_out, &samp_out);

          // random checking of the quality of our interpolations
          if (test_mode && iii%11==0 && jjj%13==0) {
            double real_line, real_samp;
            sar_to_dem(meta_sar, meta_dem, ii+iii, jj+jjj, &real_line, &real_samp);

            double err = hypot(real_line - line_out, real_samp - samp_out);

            avg_err += err;
            if (err > max_err)
              max_err = err;

            if (err > tolerance) {
              asfPrintStatus("Out of tolerance at %d,%d: (%f,%f) vs (%f,%f) -> %f\n",
                             ii+iii, jj+jjj, line_out, samp_out, real_line, real_samp,
                             err);
              ++num_out_of_tol;
            }
            if (err > .5) {
              asfPrintStatus("Error is larger than 1 pixel!\n");
              ++num_bad;
            }
            ++num_checked;
          }
          buf[index] = interp_demData(demData, dnl, dns, line_out, samp_out);
        }
      }
    }

    put_float_lines(fpOut, meta_out, ii, size, buf);
    asfPrintStatus("Completed %.1f%%  \r", 100.*ii/(double)nl);
  }
  asfPrintStatus("Completed 100%%   \n");

  if (test_mode) {
    asfPrintStatus("Tolerance was %f\n", tolerance);
    asfPrintStatus("%d/%d checked pixels had error exceeding tolerance. (%.1f%%)\n",
                   num_out_of_tol, num_checked, 100.*num_out_of_tol/(double)num_checked);
    asfPrintStatus("%d/%d checked pixels had error larger than half a pixel. (%.1f%%)\n",
                   num_bad, num_checked, 100.*num_bad/(double)num_checked);
    asfPrintStatus("Maximum error: %f pixels\n", max_err);
    avg_err /= (double)num_checked;
    asfPrintStatus("Average error: %f pixels\n", avg_err);
  }

  FCLOSE(fpOut);
  meta_write(meta_out, outImgTmp);

  meta_free(meta_out);
  meta_free(meta_dem);

  FREE(buf);
  FREE(demData);

  // now apply 3x3 filter
  if (do_averaging) {
    asfPrintStatus("Smoothing with 3x3 kernel ...\n");
    smooth(outImgTmp, outImg, 3, EDGE_TRUNCATE);
  }

  FREE(outImg);
  FREE(outImgTmp);
  FREE(output_name_tmp);

  return FALSE;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* ode-initval2/rk2imp.c
 *
 * Copyright (C) 2009, 2010 Tuomo Keskitalo
 *
 * 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.
 */

/* Based on rk2imp.c by Gerard Jungman */

/* Runge-Kutta 2, Gaussian implicit. Also known as implicit midpoint
   rule. Error estimation is carried out by the step doubling
   method.
 */

/* Reference: Ascher, U.M., Petzold, L.R., Computer methods for
   ordinary differential and differential-algebraic equations, SIAM,
   Philadelphia, 1998.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

#include "odeiv_util.h"
#include "rksubs.c"
#include "modnewton1.c"

/* Stage of method */
#define RK2IMP_STAGE 1

typedef struct
{
  gsl_matrix *A;                /* Runge-Kutta coefficients */
  double *y_onestep;            /* Result with one step */
  double *y_twostep;            /* Result with two half steps */
  double *ytmp;                 /* Temporary work space */
  double *y_save;               /* Backup space */
  double *YZ;                   /* Runge-Kutta points */
  double *fYZ;                  /* Derivatives at YZ */
  gsl_matrix *dfdy;             /* Jacobian matrix */
  double *dfdt;                 /* time derivative of f */
  modnewton1_state_t *esol;     /* nonlinear equation solver */
  double *errlev;               /* desired error level of y */
  const gsl_odeiv2_driver *driver;      /* pointer to driver object */
}
rk2imp_state_t;

static void *
rk2imp_alloc (size_t dim)
{
  rk2imp_state_t *state = (rk2imp_state_t *) malloc (sizeof (rk2imp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk2imp_state",
                      GSL_ENOMEM);
    }

  state->A = gsl_matrix_alloc (RK2IMP_STAGE, RK2IMP_STAGE);

  if (state->A == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for A", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->y_twostep = (double *) malloc (dim * sizeof (double));

  if (state->y_twostep == 0)
    {
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->y_save = (double *) malloc (dim * sizeof (double));

  if (state->y_save == 0)
    {
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_save", GSL_ENOMEM);
    }

  state->YZ = (double *) malloc (dim * RK2IMP_STAGE * sizeof (double));

  if (state->YZ == 0)
    {
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for YZ", GSL_ENOMEM);
    }

  state->fYZ = (double *) malloc (dim * RK2IMP_STAGE * sizeof (double));

  if (state->fYZ == 0)
    {
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for fYZ", GSL_ENOMEM);
    }

  state->dfdt = (double *) malloc (dim * sizeof (double));

  if (state->dfdt == 0)
    {
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdt", GSL_ENOMEM);
    }

  state->dfdy = gsl_matrix_alloc (dim, dim);

  if (state->dfdy == 0)
    {
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdy", GSL_ENOMEM);
    }

  state->esol = modnewton1_alloc (dim, RK2IMP_STAGE);

  if (state->esol == 0)
    {
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for esol", GSL_ENOMEM);
    }

  state->errlev = (double *) malloc (dim * sizeof (double));

  if (state->errlev == 0)
    {
      modnewton1_free (state->esol);
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for errlev", GSL_ENOMEM);
    }

  state->driver = NULL;

  return state;
}

static int
rk2imp_apply (void *vstate, size_t dim, double t, double h,
              double y[], double yerr[],
              const double dydt_in[], double dydt_out[],
              const gsl_odeiv2_system * sys)
{
  /* Makes a Gaussian implicit 4th order Runge-Kutta step with size h
     and estimates the local error of the step by step doubling.
   */

  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  double *const y_onestep = state->y_onestep;
  double *const y_twostep = state->y_twostep;
  double *const ytmp = state->ytmp;
  double *const y_save = state->y_save;
  double *const YZ = state->YZ;
  double *const fYZ = state->fYZ;
  gsl_matrix *const dfdy = state->dfdy;
  double *const dfdt = state->dfdt;
  double *const errlev = state->errlev;

  const modnewton1_state_t *esol = state->esol;

  /* Runge-Kutta coefficients */
  gsl_matrix *A = state->A;
  const double b[] = { 1.0 };
  const double c[] = { 0.5 };
  gsl_matrix_set (A, 0, 0, 0.5);

  if (esol == NULL)
    {
      GSL_ERROR ("no non-linear equation solver speficied", GSL_EINVAL);
    }

  /* Get desired error levels via gsl_odeiv2_control object through driver
     object, which is a requirement for this stepper.
   */

  if (state->driver == NULL)
    {
      return GSL_EFAULT;
    }
  else
    {
      size_t i;

      for (i = 0; i < dim; i++)
        {
          if (dydt_in != NULL)
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           dydt_in[i], h, i, &errlev[i]);
            }
          else
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           0.0, h, i, &errlev[i]);
            }
        }
    }

  /* Evaluate Jacobian for modnewton1 */

  {
    int s = GSL_ODEIV_JA_EVAL (sys, t, y, dfdy->data, dfdt);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Calculate a single step with size h */

  {
    int s = modnewton1_init ((void *) esol, A, h, dfdy, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h, y,
                              sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + c[0] * h, YZ, fYZ);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = rksubs (y_onestep, h, y, fYZ, b, RK2IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Error estimation by step doubling */

  {
    int s = modnewton1_init ((void *) esol, A, h / 2.0, dfdy, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* 1st half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h / 2.0, y,
                              sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      return s;
  }

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + c[0] * h / 2.0, YZ, fYZ);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = rksubs (ytmp, h / 2.0, y, fYZ, b, RK2IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Save original y values in case of error */

  DBL_MEMCPY (y_save, y, dim);

  /* 2nd half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t + h / 2.0, h / 2.0,
                              ytmp, sys, YZ, errlev);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h / 2.0 + c[0] * h / 2.0, YZ, fYZ);
    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    /* Note: rk2imp returns y using the results from two half steps
       instead of the single step since the results are freely
       available and more precise.
     */

    int s = rksubs (y_twostep, h / 2.0, ytmp, fYZ, b, RK2IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        DBL_MEMCPY (y, y_save, dim);
        return s;
      }
  }

  DBL_MEMCPY (y, y_twostep, dim);

  /* Error estimation */

  {
    size_t i;
    for (i = 0; i < dim; i++)
      {
        yerr[i] = ODEIV_ERR_SAFETY * 0.5 *
          fabs (y_twostep[i] - y_onestep[i]) / 3.0;
      }
  }

  /* Derivatives at output */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          /* Restore original values */
          DBL_MEMCPY (y, y_save, dim);

          return s;
        }
    }

  return GSL_SUCCESS;
}

static int
rk2imp_set_driver (void *vstate, const gsl_odeiv2_driver * d)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  state->driver = d;

  return GSL_SUCCESS;
}

static int
rk2imp_reset (void *vstate, size_t dim)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->y_onestep, dim);
  DBL_ZERO_MEMSET (state->y_twostep, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);
  DBL_ZERO_MEMSET (state->y_save, dim);
  DBL_ZERO_MEMSET (state->YZ, dim);
  DBL_ZERO_MEMSET (state->fYZ, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk2imp_order (void *vstate)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  state = 0;                    /* prevent warnings about unused parameters */
  return 2;
}

static void
rk2imp_free (void *vstate)
{
  rk2imp_state_t *state = (rk2imp_state_t *) vstate;

  free (state->errlev);
  modnewton1_free (state->esol);
  gsl_matrix_free (state->dfdy);
  free (state->dfdt);
  free (state->fYZ);
  free (state->YZ);
  free (state->y_save);
  free (state->ytmp);
  free (state->y_twostep);
  free (state->y_onestep);
  gsl_matrix_free (state->A);
  free (state);
}

static const gsl_odeiv2_step_type rk2imp_type = {
  "rk2imp",                     /* name */
  1,                            /* can use dydt_in? */
  1,                            /* gives exact dydt_out? */
  &rk2imp_alloc,
  &rk2imp_apply,
  &rk2imp_set_driver,
  &rk2imp_reset,
  &rk2imp_order,
  &rk2imp_free
};

const gsl_odeiv2_step_type *gsl_odeiv2_step_rk2imp = &rk2imp_type;

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef SYSTEM_H_
#define SYSTEM_H_

#pragma once

//#include <gsl/gsl_matrix.h>
#include <fstream>
//#include <gsl/gsl_vector.h>
//#include <gsl/gsl_math.h>
//#include <gsl/gsl_multimin.h>
//#include <gsl/gsl_errno.h>

#include <memory>
#include <math.h>
#include <thrust/extrema.h>
#include <thrust/sort.h>
#include <thrust/copy.h>
#include <thrust/for_each.h>
#include <thrust/transform.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/tuple.h>
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/remove.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/copy.h>
#include <thrust/execution_policy.h>
#include <thrust/pair.h>
#include <stdint.h>

//
struct CapsidInfoVecs {
	thrust::device_vector<int> id_bucket;	//bucket id
	// bucket value means global rank of a certain point
	thrust::device_vector<int> id_value;//node id

 	thrust::device_vector<double> nodeLocX;
	thrust::device_vector<double> nodeLocY;
	thrust::device_vector<double> nodeLocZ;

	//used for capside linear springs binding to membrane
	thrust::device_vector<int> tempMembraneId;
	thrust::device_vector<int> tempCapsideId;

	thrust::device_vector<double> tempLengthsPairs;
	thrust::device_vector<double> tempNodeForceX;
	thrust::device_vector<double> tempNodeForceY;
	thrust::device_vector<double> tempNodeForceZ;

	//used for capside repulsion
	int factor = 10;//number of default nodes repulsion can interact with.
	thrust::device_vector<int> tempNodeIdUnreduced;
	thrust::device_vector<double> tempNodeForceXUnreduced;
	thrust::device_vector<double> tempNodeForceYUnreduced;
	thrust::device_vector<double> tempNodeForceZUnreduced;

	thrust::device_vector<int> tempNodeIdReduced;
	thrust::device_vector<double> tempNodeForceXReduced;
	thrust::device_vector<double> tempNodeForceYReduced;
	thrust::device_vector<double> tempNodeForceZReduced;

	double spring_constant = 10.0;
	double length_zero = 0.97;
	double length_cutoff = 1.63;
	int maxNodeCount;
	double viscosity = 1.0;
	double forceX = 0.0;
	double forceY = 0.0;
	double forceZ = 0.0;

	int num_connections=0;

};

//Data Structure for node location. velocity and force
struct CoordInfoVecs {

	thrust::device_vector<double> scaling_per_edge;

	thrust::device_vector<int> nodes2Triangles_1;
	thrust::device_vector<int> nodes2Triangles_2;
	thrust::device_vector<int> nodes2Triangles_3;
	thrust::device_vector<int> nodes2Triangles_4;
	thrust::device_vector<int> nodes2Triangles_5;
	thrust::device_vector<int> nodes2Triangles_6;
	thrust::device_vector<int> nodes2Triangles_7;
	thrust::device_vector<int> nodes2Triangles_8;
	thrust::device_vector<int> nodes2Triangles_9;
//	thrust::device_vector<int> nodes2Triangles_10;
//	thrust::device_vector<int> nodes2Triangles_11;
//	thrust::device_vector<int> nodes2Triangles_12;

	thrust::device_vector<double> SurfaceNormalX;
	thrust::device_vector<double> SurfaceNormalY;
	thrust::device_vector<double> SurfaceNormalZ;

	thrust::device_vector<int> nndata1;
	thrust::device_vector<int> nndata2;
	thrust::device_vector<int> nndata3;
	thrust::device_vector<int> nndata4;
	thrust::device_vector<int> nndata5;
	thrust::device_vector<int> nndata6;
	thrust::device_vector<int> nndata7;
	thrust::device_vector<int> nndata8;
	thrust::device_vector<int> nndata9;
//	thrust::device_vector<int> nndata10;
//	thrust::device_vector<int> nndata11;
//	thrust::device_vector<int> nndata12;

	thrust::device_vector<bool> isNodeFixed;
	//GLOBAL COORDS
	// X,Y,Z, location, velocity and force of all nodes
	thrust::device_vector<double> prevNodeLocX;
	thrust::device_vector<double> prevNodeLocY;
	thrust::device_vector<double> prevNodeLocZ;

	thrust::device_vector<double> prevNodeForceX;
	thrust::device_vector<double> prevNodeForceY;
	thrust::device_vector<double> prevNodeForceZ;

 	thrust::device_vector<double> nodeLocX;
	thrust::device_vector<double> nodeLocY;
	thrust::device_vector<double> nodeLocZ;

	thrust::device_vector<double> nodeForceX;
	thrust::device_vector<double> nodeForceY;
	thrust::device_vector<double> nodeForceZ;

	thrust::device_vector<double> tempNodeForceXReduced;
	thrust::device_vector<double> tempNodeForceYReduced;
	thrust::device_vector<double> tempNodeForceZReduced;
	thrust::device_vector<double> tempNodeForceXUnreduced;
	thrust::device_vector<double> tempNodeForceYUnreduced;
	thrust::device_vector<double> tempNodeForceZUnreduced;

	//LOCAL COORDS
	//indices of each triangle
	int num_triangles;
	thrust::device_vector<int> triangles2Nodes_1;
	thrust::device_vector<int> triangles2Nodes_2;
	thrust::device_vector<int> triangles2Nodes_3;

	//indices of each edge
	int num_edges;
	thrust::device_vector<int> edges2Nodes_1;
	thrust::device_vector<int> edges2Nodes_2;

	//indices of 2 triangle on each edge
	thrust::device_vector<int> edges2Triangles_1;
	thrust::device_vector<int> edges2Triangles_2;

	//indices of edges on each triangle.
	thrust::device_vector<int> triangles2Edges_1;
	thrust::device_vector<int> triangles2Edges_2;
	thrust::device_vector<int> triangles2Edges_3;

};

//struct used for linking of nodes in network
struct AuxVecs {
	// bucket key means which bucket ID does a certain point fit into
	thrust::device_vector<int> id_bucket;	//bucket id
	// bucket value means global rank of a certain point
	thrust::device_vector<int> id_value;//node id
	// bucket key expanded means what are the bucket IDs are the neighbors of a certain point
	thrust::device_vector<int> id_bucket_expanded;
	// bucket value expanded means each point ( represented by its global rank) will have multiple copies
	thrust::device_vector<int> id_value_expanded;

	// begin position of a keys in id_bucket_expanded and id_value_expanded
	//entry keyBegin[bucketKey] returns start of indices to link
	thrust::device_vector<int> keyBegin;
	// end position of a keys in id_bucket_expanded and id_value_expanded
	thrust::device_vector<int> keyEnd;

	int endIndexid_bucket;
};

struct DomainParams {
	double minX;
	double maxX;
	double minY;
	double maxY;
	double minZ;
	double maxZ;
	double originMinX;
	double originMaxX;
	double originMinY;
	double originMaxY;
	double originMinZ;
	double originMaxZ;
	double gridSpacing = 1.5;//bucket scheme search distance, must be larger than cutoff for capsid
	int XBucketCount;
	int YBucketCount;
	int ZBucketCount;
	int totalBucketCount = 0;
};

struct LJInfoVecs{
	double LJ_PosX;
	double LJ_PosY;
	double LJ_PosZ;
	thrust::device_vector<double> LJ_PosX_all;
	thrust::device_vector<double> LJ_PosY_all;
	thrust::device_vector<double> LJ_PosZ_all;

	double Rmin_M=0.97;//1.0;
	double Rcutoff_M=0.97;//1.0;
	double Rmin_LJ;
	double Rcutoff_LJ;

	double epsilon_M=1.0;
	double epsilon_M_att1;
	double epsilon_M_att2;
	double epsilon_M_rep1;
	double epsilon_M_rep2;
	double epsilon_LJ;
	double epsilon_LJ_rep1;
	double epsilon_LJ_rep2;
	double spring_constant;

	thrust::device_vector<int> node_id_close;
	double lj_energy_M;
	double lj_energy_LJ;
	double forceX;
	double forceY;
	double forceZ;
	thrust::device_vector<double> forceX_all;
	thrust::device_vector<double> forceY_all;
	thrust::device_vector<double> forceZ_all;

};

struct AreaTriangleInfoVecs {

	int factor = 3;//used for reduction
	double initial_area = 0.0048013;
	double spring_constant;
	double spring_constant_weak;

	double area_triangle_energy;

	thrust::device_vector<int> tempNodeIdUnreduced;
	thrust::device_vector<double> tempNodeForceXUnreduced;
	thrust::device_vector<double> tempNodeForceYUnreduced;
	thrust::device_vector<double> tempNodeForceZUnreduced;

	thrust::device_vector<int> tempNodeIdReduced;
	thrust::device_vector<double> tempNodeForceXReduced;
	thrust::device_vector<double> tempNodeForceYReduced;
	thrust::device_vector<double> tempNodeForceZReduced;

};

struct BendingTriangleInfoVecs {
	int numBendingSprings=0;

	int factor = 4;//used for reduction
	double spring_constant;
	double spring_constant_weak;
	double spring_constant_raft;
	double spring_constant_coat;
	double initial_angle = 0.0;//radians
	double initial_angle_raft;
	double initial_angle_coat;

	double bending_triangle_energy;

	thrust::device_vector<int> tempNodeIdUnreduced;
	thrust::device_vector<double> tempNodeForceXUnreduced;
	thrust::device_vector<double> tempNodeForceYUnreduced;
	thrust::device_vector<double> tempNodeForceZUnreduced;

	thrust::device_vector<int> tempNodeIdReduced;
	thrust::device_vector<double> tempNodeForceXReduced;
	thrust::device_vector<double> tempNodeForceYReduced;
	thrust::device_vector<double> tempNodeForceZReduced;
};
struct LinearSpringInfoVecs {

	int factor = 2;//used for reduction
	double spring_constant;
	double spring_constant_weak;
	double spring_constant_att1;
	double spring_constant_att2;
	double spring_constant_rep1; //This is the "D" in Morse potential
	double spring_constant_rep2; //This is the "a" in Morse potential

	double linear_spring_energy;
	double memrepulsion_energy;
	double scalar_edge_length;

	thrust::device_vector<double> edge_initial_length;

	thrust::device_vector<int> tempNodeIdUnreduced;
	thrust::device_vector<double> tempNodeForceXUnreduced;
	thrust::device_vector<double> tempNodeForceYUnreduced;
	thrust::device_vector<double> tempNodeForceZUnreduced;

	thrust::device_vector<int> tempNodeIdReduced;
	thrust::device_vector<double> tempNodeForceXReduced;
	thrust::device_vector<double> tempNodeForceYReduced;
	thrust::device_vector<double> tempNodeForceZReduced;
};

struct GeneralParams{
	double kT;
	double kT_growth;
	double tau;
	int solve_time=100;
	double Rmin = 1.0; //Current mesh minimum edge length, subject to change.
	double abs_Rmin;
	int iteration = 0;
	int maxNodeCount;
	int maxNodeCountLJ;
	double length_scale;
	//parameters for advancing timestep and determining equilibrium

	double dt;
	double nodeMass = 1.0;
	int edges_in_upperhem_list_length;
	double growth_energy_scaling;
	thrust::device_vector<int> edge_to_ljparticle;
	thrust::device_vector<int> nodes_in_upperhem;
	thrust::device_vector<int> edges_in_upperhem;
	thrust::device_vector<int> edges_in_upperhem_list;
	thrust::device_vector<int> triangles_in_upperhem;
	thrust::device_vector<int> boundaries_in_upperhem;
	thrust::device_vector<double> angle_per_edge;
	double centerX = 0.0;
	double centerY = 0.0;
	double centerZ = 0.0;

	double current_total_volume;
	double true_current_total_volume;
	double eq_total_volume;
	double volume_spring_constant;
	double volume_energy;
	double eq_total_boundary_length;
	double line_tension_energy;
	double line_tension_constant;
	double safeguardthreshold;

	//int num_of_triangles;
	//int num_of_edges;
	int true_num_edges;

	double insertion_energy_cost;
	double strain_threshold;

	int SCALE_TYPE;
	double scaling_pow;
	double gausssigma;
	double hilleqnconst;
	double hilleqnpow;

	thrust::device_vector<int> no_weakening;
	double septin_ring_z;
	double boundary_z;

};

class Storage;
class SystemBuilder;
struct HostSetInfoVecs;

class System {
public:
	std::weak_ptr<SystemBuilder> weak_bld_ptr;
	GeneralParams generalParams;
	DomainParams domainParams;
	AuxVecs auxVecs;
	CoordInfoVecs coordInfoVecs;

	CapsidInfoVecs capsidInfoVecs;
	LinearSpringInfoVecs linearSpringInfoVecs;
	BendingTriangleInfoVecs bendingTriangleInfoVecs;
	AreaTriangleInfoVecs areaTriangleInfoVecs;
	LJInfoVecs ljInfoVecs;

	std::shared_ptr<Storage> storage;

	//gsl_vector* df;
	//gsl_vector* locations;

public:

	System();

	void set_weak_builder(std::weak_ptr<SystemBuilder> _weak_bld_ptr);

	void PrintForce();

	void initializeSystem(HostSetInfoVecs& hostSetInfoVecs);

	void assignStorage(std::shared_ptr<Storage> _storage);

	void solveSystem();

	void setExtras();

	void setBucketScheme();

	//void Solve_Forces(const gsl_vector* temp_locations);
	void Solve_Forces();
	//double Solve_Energy(const gsl_vector* temp_locations);

	//void dev_to_gsl_loc_update(gsl_vector* temp_locations);
	//void gsl_to_dev_loc_update(const gsl_vector* temp_locations);

};

#endif /*POLYMERSYSTEM_H_*/


text:
#ifndef __GSL_PERMUTE_MATRIX_H__
#define __GSL_PERMUTE_MATRIX_H__

#include <gsl/gsl_permute_matrix_complex_long_double.h>
#include <gsl/gsl_permute_matrix_complex_double.h>
#include <gsl/gsl_permute_matrix_complex_float.h>

#include <gsl/gsl_permute_matrix_long_double.h>
#include <gsl/gsl_permute_matrix_double.h>
#include <gsl/gsl_permute_matrix_float.h>

#include <gsl/gsl_permute_matrix_ulong.h>
#include <gsl/gsl_permute_matrix_long.h>

#include <gsl/gsl_permute_matrix_uint.h>
#include <gsl/gsl_permute_matrix_int.h>

#include <gsl/gsl_permute_matrix_ushort.h>
#include <gsl/gsl_permute_matrix_short.h>

#include <gsl/gsl_permute_matrix_uchar.h>
#include <gsl/gsl_permute_matrix_char.h>

#endif /* __GSL_PERMUTE_MATRIX_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest


text:
#ifndef DEPTH_TO_BEARING_H_ZXFA9HGG
#define DEPTH_TO_BEARING_H_ZXFA9HGG

#include <cmath>
#include <gsl/gsl>
#include <limits>
#include <sens_loc/camera_models/concepts.h>
#include <sens_loc/camera_models/utility.h>
#include <sens_loc/conversion/util.h>
#include <sens_loc/math/constants.h>
#include <sens_loc/math/coordinate.h>
#include <sens_loc/math/image.h>
#include <sens_loc/math/triangles.h>
#include <sens_loc/util/correctness_util.h>
#include <taskflow/taskflow.hpp>

namespace sens_loc::conversion {

/// Convert the image \p depth_image to an bearing angle image.
/// This function returns a new image with the same dimension as \p depth_image.
///
/// \tparam Direction the direction for neighbourhood-relationship between
/// pixels that form the bearing angle
/// \tparam Real precision to calculate in, needs to be floating-point
/// \tparam Intrinsic camera model that projects pixel to the unit sphere
/// \param depth_image range image that was taken by a sensor with the
/// calibration from \p intrinsic
/// \param intrinsic camera model to calculate the angle between light
/// rays that correspond to pixels
/// \returns cv::Mat<Real> with each pixel the bearing angle (radians) in
/// \p Direction. Invalid depth values result in 0 as result.
//
/// \note Depth images are orthografic and require conversion first!
/// \sa conversion::depth_to_laserscan
/// \sa camera_models::is_intrinsic_v
///
/// \pre \p depth_image to have 1 channel
/// \pre \p depth_image == laser-scan like image!
/// \pre \p depth_image is not empty
/// \pre \p the underlying type of \p depth_image matches \p PixelType
/// \pre \p intrinsic matches sensor that took the image, prior preprocessing
/// \post each bearing angle is in the range \f$[0, \pi)\f$
/// is of course possible with matching changes in the \p depth_image.
template <direction Direction,
          template <typename>
          typename Intrinsic,
          typename Real = float>
math::image<Real> depth_to_bearing(const math::image<Real>& depth_image,
                                   const Intrinsic<Real>&   intrinsic) noexcept;

/// This function provides are parallelized version of the conversion
/// functionality.
///
/// This function creates a taskflow for the row-wise parallel calculation
/// of the bearing angle image.
/// Only differences are documented here.
/// \param[in] depth_image,intrinsic the same
/// \param[out] ba_image result image that will be created with parallel
/// processing
/// \param[inout] flow parallel flow type that is used to parallelize the outer
/// for loop over all rows.
/// \returns synchronization points before and after the calculation of the
/// bearing angle image.
/// \sa depth_to_bearing
template <direction Direction,
          template <typename>
          typename Intrinsic,
          typename Real = float>
std::pair<tf::Task, tf::Task>
par_depth_to_bearing(const math::image<Real>& depth_image,
                     const Intrinsic<Real>&   intrinsic,
                     math::image<Real>&       ba_image,
                     tf::Taskflow&            flow) noexcept;

/// Convert a bearing angle image to an image with integer types.
/// This function scales the bearing angles between
/// [PixelType::min, PixelType::max] for the angles in range (0, PI).
///
/// \tparam Real underyling type of the \p bearing_image, floating-point
/// \tparam PixelType target type for the converted image that is returned,
/// arithmetic type (integer or floating-point)
/// \param bearing_image calculated
/// bearing angle image
/// \returns the bearing angle image is converted to a "normal" image that can
/// be displayed and stored with normal image visualization tools.
/// \pre the underlying type of \p bearing_image is \p Real
/// \pre range of pixel is in \p bearinge_image \f$[0, \pi)\f$
/// \post the underlying type of the result is \p PixelType
/// \post range of pixel in result is \f$[PixelType_{min}, PixelType_{max}]\f$
/// \sa conversion::depth_to_bearing
template <typename Real = float, typename PixelType = ushort>
math::image<PixelType>
convert_bearing(const math::image<Real>& bearing_image) noexcept;

namespace detail {
inline int get_du(direction dir) {
    switch (dir) {
    case direction::antidiagonal:
    case direction::diagonal:
    case direction::horizontal: return -1;
    case direction::vertical: return 0;
    }
    UNREACHABLE("Only 4 directions are possible");  // LCOV_EXCL_LINE
}

inline int get_dv(direction dir) {
    switch (dir) {
    case direction::diagonal:
    case direction::vertical: return -1;
    case direction::horizontal: return 0;
    case direction::antidiagonal: return 1;
    }
    UNREACHABLE("Only 4 directions are possible");  // LCOV_EXCL_LINE
}

/// Provide a general accessor for the pixels that provide depth information
/// to calculate the bearing angle.
template <typename Real, direction Direction>
class pixel {
  public:
    pixel()
        : _du{get_du(Direction)}
        , _dv{get_dv(Direction)} {}
    /// \returns (u, v) for the prior point depending on the direction.
    math::pixel_coord<int> operator()(const math::pixel_coord<int>& p) const
        noexcept {
        return {p.u() + _du, p.v() + _dv};
    }

  private:
    const int _du = 0;
    const int _dv = 0;
};

inline int get_x_start(direction dir) {
    switch (dir) {
    case direction::horizontal:
    case direction::antidiagonal:
    case direction::diagonal: return 1;
    case direction::vertical: return 0;
    }
    UNREACHABLE("Only 4 directions are possible");  // LCOV_EXCL_LINE
}

inline int get_y_start(direction dir) {
    switch (dir) {
    case direction::horizontal:
    case direction::antidiagonal: return 0;
    case direction::diagonal:
    case direction::vertical: return 1;
    }
    UNREACHABLE("Only 4 directions are possible");  // LCOV_EXCL_LINE
}

inline int get_dy_end(direction dir) {
    switch (dir) {
    case direction::antidiagonal: return 1;
    case direction::horizontal:
    case direction::diagonal:
    case direction::vertical: return 0;
    }
    UNREACHABLE("Only 4 directions are possible");  // LCOV_EXCL_LINE
}

/// Define the start and end iterator for each bearing angle direction.
template <direction Direction>
struct pixel_range {
    pixel_range(const cv::Mat& depth_image) noexcept
        : x_start{get_x_start(Direction)}
        , x_end{depth_image.cols}
        , y_start{get_y_start(Direction)}
        , y_end{depth_image.rows - get_dy_end(Direction)} {}

    const int x_start;
    const int x_end;
    const int y_start;
    const int y_end;
};

template <typename Real,
          typename RangeLimits,
          typename PriorAccess,
          template <typename>
          typename Intrinsic>
inline void bearing_inner(const RangeLimits&       r,
                          const PriorAccess&       prior_accessor,
                          const int                v,
                          const math::image<Real>& depth_image,
                          const Intrinsic<Real>&   intrinsic,
                          math::image<Real>&       ba_image) {
    for (int u = r.x_start; u < r.x_end; ++u) {
        const math::pixel_coord<int> central(u, v);
        const math::pixel_coord<int> prior = prior_accessor(central);

        const auto d_i = depth_image.at(central);
        const auto d_j = depth_image.at(prior);

        Expects(d_i >= Real(0.));
        Expects(d_j >= Real(0.));

        const Real d_phi = camera_models::phi(intrinsic, central, prior);

        // A depth==0 means there is no measurement at this pixel.
        const Real angle =
            (d_i == Real(0.) || d_j == Real(0.))
                ? Real(0.)
                : math::bearing_angle<Real>(d_i, d_j, std::cos(d_phi));

        Ensures(angle >= Real(0.));
        Ensures(angle < math::pi<Real>);

        ba_image.at(central) = angle;
    }
}
}  // namespace detail

template <direction Direction,
          template <typename>
          typename Intrinsic,
          typename Real>
inline math::image<Real>
depth_to_bearing(const math::image<Real>& depth_image,
                 const Intrinsic<Real>&   intrinsic) noexcept {
    static_assert(camera_models::is_intrinsic_v<Intrinsic, Real>);
    static_assert(std::is_floating_point_v<Real>);

    Expects(depth_image.w() == intrinsic.w());
    Expects(depth_image.h() == intrinsic.h());

    using namespace detail;
    const pixel<Real, Direction> prior_accessor;
    const pixel_range<Direction> r{depth_image.data()};

    // Image of Reals, that will be converted after the full calculation.
    cv::Mat ba(depth_image.h(), depth_image.w(),
               math::detail::get_opencv_type<Real>());
    ba = Real(0.);
    math::image<Real> ba_image(std::move(ba));

    for (int v = r.y_start; v < r.y_end; ++v)
        detail::bearing_inner(r, prior_accessor, v, depth_image, intrinsic,
                              ba_image);

    Ensures(ba_image.h() == depth_image.h());
    Ensures(ba_image.w() == depth_image.w());

    return ba_image;
}

template <direction Direction,
          template <typename>
          typename Intrinsic,
          typename Real>
inline std::pair<tf::Task, tf::Task>
par_depth_to_bearing(const math::image<Real>& depth_image,
                     const Intrinsic<Real>&   intrinsic,
                     math::image<Real>&       ba_image,
                     tf::Taskflow&            flow) noexcept {
    static_assert(camera_models::is_intrinsic_v<Intrinsic, Real>);
    static_assert(std::is_floating_point_v<Real>);

    using namespace detail;
    Expects(depth_image.w() == intrinsic.w());
    Expects(depth_image.h() == intrinsic.h());
    Expects(ba_image.w() == depth_image.w());
    Expects(ba_image.h() == depth_image.h());

    const pixel<Real, Direction> prior_accessor;
    const pixel_range<Direction> r{depth_image.data()};

    auto sync_points = flow.parallel_for(
        r.y_start, r.y_end, 1,
        [prior_accessor, r, &depth_image, &intrinsic, &ba_image](int v) {
            detail::bearing_inner(r, prior_accessor, v, depth_image, intrinsic,
                                  ba_image);
        });

    return sync_points;
}

template <typename Real, typename PixelType>
inline math::image<PixelType>
convert_bearing(const math::image<Real>& bearing_image) noexcept {
    static_assert(std::is_floating_point_v<Real>);
    static_assert(std::is_arithmetic_v<PixelType>);

    using detail::scaling_factor;
    cv::Mat img(bearing_image.h(), bearing_image.w(),
                math::detail::get_opencv_type<PixelType>());
    auto [scale, offset] =
        scaling_factor<Real, PixelType>(/*max_angle = */ math::pi<Real>);
    bearing_image.data().convertTo(
        img, math::detail::get_opencv_type<PixelType>(), scale, offset);

    Ensures(img.cols == bearing_image.w());
    Ensures(img.rows == bearing_image.h());
    Ensures(img.type() == math::detail::get_opencv_type<PixelType>());
    Ensures(img.channels() == 1);

    return math::image<PixelType>(std::move(img));
}
}  // namespace sens_loc::conversion

#endif /* end of include guard: DEPTH_TO_BEARING_H_ZXFA9HGG */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions