text:
#ifndef NOVUGNO_IMI_INCLUDE_GNOIMI_H_
#define NOVUGNO_IMI_INCLUDE_GNOIMI_H_

#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <thread>
#include <functional>
#include <vector>
#include <ctime>
#include <chrono>

#include <mutex>
#include <cblas.h>

struct threadInfo;

class GNOIMI {
    public:
    private:
	int D;
	int K;
	int totalLearnCount;
	int learnIterationsCount;
	int L;

	int threadsCount;
	int totalPoints_perThread;
	int Init_chunkSize;

	int* coarseAssigns;
	int* fineAssigns;
	float* alphaNum;
	float* alphaDen;
	float* alpha;

	vector<float*> alphaNumerators;
	vector<float*> alphaDenominators;

	float* coarseVocab;
	float* fineVocab;
	float* fineVocabNum;
	float* fineVocabDen;
	float* coarseVocabNum;
	float* coarseVocabDen;

	vector<float*> fineVocabNumerators;
	vector<float*> fineVocabDenominators;
	vector<float*> coarseVocabNumerators;
	vector<float*> coarseVocabDenominators;

	float* coarseNorms;
	float* fineNorms;
	float* coarseFineProducts;

	float* errors;

	string output_prefix;
};

using namespace std;

int D = 512;
int K = 657;
int totalLearnCount = 260000;
int learnIterationsCount = 20;
int L = 6;

string learnFilename = "./normalized_features_260k.fvecs";
string initCoarseFilename = "./coarse_260k.fvecs";
string initFineFilename = "./fine_260k.fvecs";
string outputFilesPrefix = "./260K_new_";

int threadsCount = 3;
int totalPoints_perThread = 0;
const int Init_chunkSize = 3000;

int* coarseAssigns = (int*)malloc(totalLearnCount * sizeof(int));
int* fineAssigns = (int*)malloc(totalLearnCount * sizeof(int));
float* alphaNum = (float*)malloc(K * K * sizeof(float));
float* alphaDen = (float*)malloc(K * K * sizeof(float));
float* alpha = (float*)malloc(K * K * sizeof(float));

vector<float*> alphaNumerators(threadsCount);
vector<float*> alphaDenominators(threadsCount);

float* coarseVocab = (float*)malloc(D * K * sizeof(float));
float* fineVocab = (float*)malloc(D * K * sizeof(float));
float* fineVocabNum = (float*)malloc(D * K * sizeof(float));
float* fineVocabDen = (float*)malloc(K * sizeof(float));
float* coarseVocabNum = (float*)malloc(D * K * sizeof(float));
float* coarseVocabDen = (float*)malloc(K * sizeof(float));

vector<float*> fineVocabNumerators(threadsCount);
vector<float*> fineVocabDenominators(threadsCount);
vector<float*> coarseVocabNumerators(threadsCount);
vector<float*> coarseVocabDenominators(threadsCount);

float* coarseNorms = (float*)malloc(K * sizeof(float));
float* fineNorms = (float*)malloc(K * sizeof(float));
float* coarseFineProducts = (float*)malloc(K * K * sizeof(float));

float* errors = (float*)malloc(threadsCount * sizeof(float));
mutex mtx;

struct threadInfo {
    long long startId;
    int pointsCount;
    int chunksCount;
    int trainThreadChunkSize;
    int extra_ChunkSize;
};

threadInfo get_threadInfo(const int& threadId) {
    threadInfo t_info;
    t_info.startId = totalPoints_perThread * threadId;
    t_info.trainThreadChunkSize = Init_chunkSize;

    if (totalPoints_perThread * (threadId + 1) <= totalLearnCount) {
	t_info.pointsCount = totalPoints_perThread;
    } else {
	t_info.pointsCount = totalLearnCount - totalPoints_perThread * threadId;
    }

    t_info.extra_ChunkSize = 0;
    if (t_info.pointsCount % t_info.trainThreadChunkSize == 0) {
	t_info.chunksCount = t_info.pointsCount / t_info.trainThreadChunkSize;
    } else {
	t_info.extra_ChunkSize = t_info.pointsCount % t_info.trainThreadChunkSize;
	t_info.chunksCount = (t_info.pointsCount - t_info.extra_ChunkSize) / t_info.trainThreadChunkSize + 1;
    }

    return t_info;
}

void computeOptimalAssignsSubset(int threadId) {

    threadInfo cur_t_info = get_threadInfo(threadId);
    int trainThreadChunkSize = cur_t_info.trainThreadChunkSize;

    errors[threadId] = 0.0;

    FILE* learnStream = fopen(learnFilename.c_str(), "r");
    fseek(learnStream, cur_t_info.startId * (D + 1) * sizeof(float), SEEK_SET);

    float* pointsCoarseTerms = (float*)malloc(trainThreadChunkSize * K * sizeof(float));
    float* pointsFineTerms = (float*)malloc(trainThreadChunkSize * K * sizeof(float));
    float* chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));

    std::vector<std::pair<float, int> > coarseScores(K);

    for(int chunkId = 0; chunkId < cur_t_info.chunksCount; ++chunkId) {
	std::cout << "[Assigns][Thread " << threadId << "] " << "processing chunk " <<  chunkId << " of " << cur_t_info.chunksCount << "\n";

	if (chunkId == cur_t_info.chunksCount - 1 && cur_t_info.extra_ChunkSize != 0) {
	    free(pointsCoarseTerms);
	    free(pointsFineTerms);
	    free(chunkPoints);

	    trainThreadChunkSize = cur_t_info.extra_ChunkSize;
	    pointsCoarseTerms = (float*)malloc(trainThreadChunkSize * K * sizeof(float));
	    pointsFineTerms = (float*)malloc(trainThreadChunkSize * K * sizeof(float));
	    chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));
	}

	fvecs_fread(learnStream, chunkPoints, trainThreadChunkSize, D);
	fmat_mul_full(coarseVocab, chunkPoints, K, trainThreadChunkSize, D, "TN", pointsCoarseTerms);
	fmat_mul_full(fineVocab, chunkPoints, K, trainThreadChunkSize, D, "TN", pointsFineTerms);

	for(int pointId = 0; pointId < trainThreadChunkSize; ++pointId) {
	    cblas_saxpy(K, -1.0, coarseNorms, 1, pointsCoarseTerms + pointId * K, 1);

	    for(int k = 0; k < K; ++k) {
		coarseScores[k].first = (-1.0) * pointsCoarseTerms[pointId * K + k];
		coarseScores[k].second = k;
	    }

	    std::sort(coarseScores.begin(), coarseScores.end());
	    float currentMinScore = 999999999.0;
	    int currentMinCoarseId = -1;
	    int currentMinFineId = -1;

	    for(int l = 0; l < L; ++l) {
		//examine cluster l
		int currentCoarseId = coarseScores[l].second;
		float currentCoarseTerm = coarseScores[l].first;
		for(int currentFineId = 0; currentFineId < K; ++currentFineId) {
		    float alphaFactor = alpha[currentCoarseId * K + currentFineId];
		    float score = currentCoarseTerm + alphaFactor * coarseFineProducts[currentCoarseId * K + currentFineId] +
			(-1.0) * alphaFactor * pointsFineTerms[pointId * K + currentFineId] +
			alphaFactor * alphaFactor * fineNorms[currentFineId];

		    if(score < currentMinScore) {
			currentMinScore = score;
			currentMinCoarseId = currentCoarseId;
			currentMinFineId = currentFineId;
		    }
		}
	    }

	    coarseAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId] = currentMinCoarseId;
	    fineAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId] = currentMinFineId;

	    errors[threadId] += currentMinScore * 2 + 1.0; // point has a norm equals 1.0
	}
    }

    fclose(learnStream);
    free(chunkPoints);
    free(pointsCoarseTerms);
    free(pointsFineTerms);
}

void computeOptimalAlphaSubset(int threadId) {
    memset(alphaNumerators[threadId], 0, K * K * sizeof(float));
    memset(alphaDenominators[threadId], 0, K * K * sizeof(float));

    threadInfo cur_t_info = get_threadInfo(threadId);
    int trainThreadChunkSize = cur_t_info.trainThreadChunkSize;

    FILE* learnStream = fopen(learnFilename.c_str(), "r");
    fseek(learnStream, cur_t_info.startId * (D + 1) * sizeof(float), SEEK_SET);

    float* residual = (float*)malloc(D * sizeof(float));
    float* chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));

    for(int chunkId = 0; chunkId < cur_t_info.chunksCount; ++chunkId) {
	std::cout << "[Alpha][Thread " << threadId << "] " << "processing chunk " <<  chunkId << " of " << cur_t_info.chunksCount << "\n";

	if (chunkId == cur_t_info.chunksCount - 1 && cur_t_info.extra_ChunkSize != 0) {
	    free(chunkPoints);

	    trainThreadChunkSize = cur_t_info.extra_ChunkSize;
	    chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));
	}

	fvecs_fread(learnStream, chunkPoints, trainThreadChunkSize, D);
	for(int pointId = 0; pointId < trainThreadChunkSize; ++pointId) {
	    //S_k
	    int coarseAssign = coarseAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];
	    //T_l
	    int fineAssign = fineAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];

	    memcpy(residual, chunkPoints + pointId * D, D * sizeof(float));
	    //P_i - S_k
	    cblas_saxpy(D, -1.0, coarseVocab + coarseAssign * D, 1, residual, 1);

	    //<P_i - S_k, T_l>
	    alphaNumerators[threadId][coarseAssign * K + fineAssign] +=
		cblas_sdot(D, residual, 1, fineVocab + fineAssign * D, 1);

	    //||T||^2
	    alphaDenominators[threadId][coarseAssign * K + fineAssign] += fineNorms[fineAssign] * 2; // we keep halves of norms
	}
    }
    fclose(learnStream);
    free(chunkPoints);
    free(residual);
}

void computeOptimalFineVocabSubset(int threadId) {
    memset(fineVocabNumerators[threadId], 0, K * D * sizeof(float));
    memset(fineVocabDenominators[threadId], 0, K * sizeof(float));

    threadInfo cur_t_info = get_threadInfo(threadId);
    int trainThreadChunkSize = cur_t_info.trainThreadChunkSize;

    FILE* learnStream = fopen(learnFilename.c_str(), "r");
    fseek(learnStream, cur_t_info.startId * (D + 1) * sizeof(float), SEEK_SET);

    float* residual = (float*)malloc(D * sizeof(float));
    float* chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));

    for(int chunkId = 0; chunkId < cur_t_info.chunksCount; ++chunkId) {
	//std::cout << "[Fine vocabs][Thread " << threadId << "] " << "processing chunk " <<  chunkId << " of " << chunksCount << "\n";

	if (chunkId == cur_t_info.chunksCount - 1 && cur_t_info.extra_ChunkSize != 0) {
	    free(chunkPoints);

	    trainThreadChunkSize = cur_t_info.extra_ChunkSize;
	    chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));
	}

	fvecs_fread(learnStream, chunkPoints, trainThreadChunkSize, D);
	for(int pointId = 0; pointId < trainThreadChunkSize; ++pointId) {

	    int coarseAssign = coarseAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];
	    int fineAssign = fineAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];

	    float alphaFactor = alpha[coarseAssign * K + fineAssign];
	    memcpy(residual, chunkPoints + pointId * D, D * sizeof(float));

	    cblas_saxpy(D, -1.0, coarseVocab + coarseAssign * D, 1, residual, 1);
	    cblas_saxpy(D, alphaFactor, residual, 1, fineVocabNumerators[threadId] + fineAssign * D, 1);
	    fineVocabDenominators[threadId][fineAssign] += alphaFactor * alphaFactor;
	}
    }

    fclose(learnStream);
    free(chunkPoints);
    free(residual);
}

void computeOptimalCoarseVocabSubset(int threadId) {
    memset(coarseVocabNumerators[threadId], 0, K * D * sizeof(float));
    memset(coarseVocabDenominators[threadId], 0, K * sizeof(float));

    threadInfo cur_t_info = get_threadInfo(threadId);
    int trainThreadChunkSize = cur_t_info.trainThreadChunkSize;

    FILE* learnStream = fopen(learnFilename.c_str(), "r");
    fseek(learnStream, cur_t_info.startId * (D + 1) * sizeof(float), SEEK_SET);

    float* residual = (float*)malloc(D * sizeof(float));
    float* chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));

    for(int chunkId = 0; chunkId < cur_t_info.chunksCount; ++chunkId) {
	//std::cout << "[Coarse vocabs][Thread " << threadId << "] " << "processing chunk " <<  chunkId << " of " << chunksCount << "\n";

	if (chunkId == cur_t_info.chunksCount - 1 && cur_t_info.extra_ChunkSize != 0) {
	    free(chunkPoints);

	    trainThreadChunkSize = cur_t_info.extra_ChunkSize;
	    chunkPoints = (float*)malloc(trainThreadChunkSize * D * sizeof(float));
	}

	fvecs_fread(learnStream, chunkPoints, trainThreadChunkSize, D);
	for(int pointId = 0; pointId < trainThreadChunkSize; ++pointId) {
	    int coarseAssign = coarseAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];
	    int fineAssign = fineAssigns[cur_t_info.startId + chunkId * cur_t_info.trainThreadChunkSize + pointId];
	    float alphaFactor = alpha[coarseAssign * K + fineAssign];
	    memcpy(residual, chunkPoints + pointId * D, D * sizeof(float));
	    cblas_saxpy(D, -1.0 * alphaFactor, fineVocab + fineAssign * D, 1, residual, 1);
	    cblas_saxpy(D, 1, residual, 1, coarseVocabNumerators[threadId] + coarseAssign * D, 1);
	    coarseVocabDenominators[threadId][coarseAssign] += 1.0;
	}
    }

    fclose(learnStream);
    free(chunkPoints);
    free(residual);
}

int main() {
    totalPoints_perThread = totalLearnCount / threadsCount;

    //increase the threadCount
    if (totalLearnCount % threadsCount != 0) {
       	threadsCount += 1;
	alphaNumerators.resize(threadsCount);
	alphaDenominators.resize(threadsCount);
	fineVocabNumerators.resize(threadsCount);
	fineVocabDenominators.resize(threadsCount);
	coarseVocabNumerators.resize(threadsCount);
	coarseVocabDenominators.resize(threadsCount);
	free(errors);
	errors = (float*)malloc(threadsCount * sizeof(float));
    }

    for(int threadId = 0; threadId < threadsCount; ++threadId) {
	alphaNumerators[threadId] = (float*)malloc(K * K * sizeof(float*));
	alphaDenominators[threadId] = (float*)malloc(K * K * sizeof(float*));
    }

    for(int threadId = 0; threadId < threadsCount; ++threadId) {
	fineVocabNumerators[threadId] = (float*)malloc(K * D * sizeof(float*));
	fineVocabDenominators[threadId] = (float*)malloc(K * sizeof(float*));
    }

    for(int threadId = 0; threadId < threadsCount; ++threadId) {
	coarseVocabNumerators[threadId] = (float*)malloc(K * D * sizeof(float));
	coarseVocabDenominators[threadId] = (float*)malloc(K * sizeof(float));
    }

    // init vocabs
    fvecs_read(initCoarseFilename.c_str(), D, K, coarseVocab);
    fvecs_read(initFineFilename.c_str(), D, K, fineVocab);

    // init alpha
    for(int i = 0; i < K * K; ++i) {
	alpha[i] = 1.0f;
    }

    // learn iterations
    std::cout << "Start learning iterations...\n";
    for(int it = 0; it < learnIterationsCount; ++it) {
	cout << "learning iteration " << it << endl;
	for(int k = 0; k < K; ++k) {
	    coarseNorms[k] = cblas_sdot(D, coarseVocab + k * D, 1, coarseVocab + k * D, 1) / 2;
	    fineNorms[k] = cblas_sdot(D, fineVocab + k * D, 1, fineVocab + k * D, 1) / 2;
	}

	fmat_mul_full(fineVocab, coarseVocab, K, K, D, "TN", coarseFineProducts);

	// update Assigns
	vector<std::thread> workers;
	memset(errors, 0, threadsCount * sizeof(float));
	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers.push_back(std::thread(computeOptimalAssignsSubset, threadId));
	}

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers[threadId].join();
	}

	float totalError = 0.0;
	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    cout << "threadId " << threadId << ": err " << errors[threadId] << endl;
	    totalError += errors[threadId];
	}
	std::cout << "Current reconstruction error... " << totalError / totalLearnCount << "\n";

	// update alpha
	workers.clear();
	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers.push_back(std::thread(computeOptimalAlphaSubset, threadId));
	}
	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers[threadId].join();
	}

	// update fine Vocabs
	workers.clear();
	memset(alphaNum, 0, K * K * sizeof(float));
	memset(alphaDen, 0, K * K * sizeof(float));

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    cblas_saxpy(K * K, 1, alphaNumerators[threadId], 1, alphaNum, 1);
	    cblas_saxpy(K * K, 1, alphaDenominators[threadId], 1, alphaDen, 1);
	}

	for(int i = 0; i < K * K; ++i) {
	    alpha[i] = (alphaDen[i] == 0) ? 1.0 : alphaNum[i] / alphaDen[i];
	}

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers.push_back(std::thread(computeOptimalFineVocabSubset, threadId));
	}

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers[threadId].join();
	}

	// update coarse Vocabs
	workers.clear();
	memset(fineVocabNum, 0, K * D * sizeof(float));
	memset(fineVocabDen, 0, K * sizeof(float));

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    cblas_saxpy(K * D, 1, fineVocabNumerators[threadId], 1, fineVocabNum, 1);
	    cblas_saxpy(K, 1, fineVocabDenominators[threadId], 1, fineVocabDen, 1);
	}

	for(int i = 0; i < K * D; ++i) {
	    fineVocab[i] = (fineVocabDen[i / D] == 0) ? 0 : fineVocabNum[i] / fineVocabDen[i / D];
	}

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers.push_back(std::thread(computeOptimalCoarseVocabSubset, threadId));
	}

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    workers[threadId].join();
	}

	workers.clear();
	memset(coarseVocabNum, 0, K * D * sizeof(float));
	memset(coarseVocabDen, 0, K * sizeof(float));

	for(int threadId = 0; threadId < threadsCount; ++threadId) {
	    cblas_saxpy(K * D, 1, coarseVocabNumerators[threadId], 1, coarseVocabNum, 1);
	    cblas_saxpy(K, 1, coarseVocabDenominators[threadId], 1, coarseVocabDen, 1);
	}

	for(int i = 0; i < K * D; ++i) {
	    coarseVocab[i] = (coarseVocabDen[i / D] == 0) ? 0 : coarseVocabNum[i] / coarseVocabDen[i / D];
	}

	// save current alpha and vocabs
	std::stringstream alphaFilename;
	alphaFilename << outputFilesPrefix << "alpha_" << it << ".fvecs";
	fvecs_write(alphaFilename.str().c_str(), K, K, alpha);

	std::stringstream fineVocabFilename;
	fineVocabFilename << outputFilesPrefix << "fine_" << it << ".fvecs";
	fvecs_write(fineVocabFilename.str().c_str(), D, K, fineVocab);

	std::stringstream coarseVocabFilename;
	coarseVocabFilename << outputFilesPrefix << "coarse_" << it << ".fvecs";
	fvecs_write(coarseVocabFilename.str().c_str(), D, K, coarseVocab);

	cout << "iteration " << it << " finished!" << endl;
    }

    free(coarseAssigns);
    free(fineAssigns);
    free(alphaNum);
    free(alphaDen);
    free(alpha);
    free(coarseVocab);
    free(coarseVocabNum);
    free(coarseVocabDen);
    free(fineVocab);
    free(fineVocabNum);
    free(fineVocabDen);
    free(coarseNorms);
    free(fineNorms);
    free(coarseFineProducts);
    free(errors);
    return 0;
}
#endif  //NOVUGNO_IMI_INCLUDE_GNOIMI_H_


text:
#ifndef KEYPOINT_TRANSFORM_H_ARXFM8VJ
#define KEYPOINT_TRANSFORM_H_ARXFM8VJ

#include <gsl/gsl>
#include <sens_loc/camera_models/pinhole.h>
#include <sens_loc/camera_models/projection.h>
#include <sens_loc/math/pointcloud.h>

namespace sens_loc::apps {

template <template <typename> typename Model = sens_loc::camera_models::pinhole,
          typename Real                      = float>
math::pointcloud<Real>
keypoints_to_pointcloud(const std::vector<cv::KeyPoint>& kps,
                        const math::image<ushort>&       depth_image,
                        const Model<Real>&               c,
                        float                            unit_factor) noexcept {
    using camera_models::keypoint_to_coords;
    math::imagepoints<Real> pts   = keypoint_to_coords(kps);
    math::pointcloud<Real>  c_pts = project_to_sphere(c, pts);

    for (std::size_t i = 0; i < c_pts.size(); ++i) {
        auto orig_depth = depth_image.at(pts[i]);
        auto t_d        = unit_factor * gsl::narrow_cast<float>(orig_depth);
        c_pts[i]        = t_d * c_pts[i];
    }

    return c_pts;
}

/// Transform the pointcloud \c points_in_other_frame into the coordinate
/// frame of \c c and project these points into \c c's image plane.
/// \returns the backprojected points
/// \sa project_to_image
template <template <typename> typename Model = sens_loc::camera_models::pinhole,
          typename Real                      = float>
math::imagepoints<Real>
project_to_other_camera(const math::pose_t&           pose,
                        const math::pointcloud<Real>& points_in_other_frame,
                        const Model<Real>&            c) noexcept {
    math::pointcloud<Real>  transformed   = pose * points_in_other_frame;
    math::imagepoints<Real> backprojected = project_to_image(c, transformed);
    return backprojected;
}

}  // namespace sens_loc::apps

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

ext: h

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
/*
** Copyright (C) 2004 Jonathan G. Underwood <j.underwood@open.ac.uk>
**
** 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.
*/

/* Things to check: Not sure I'm using the correct error macros - should I be
   using the ones in the GSL manual? Based the ones here on what i see in
   coupling.c. */

#include <config.h>
#include <stdlib.h>
#include <error.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_wigner.h>

/* Static prototypes. */
static int istriangle (const int two_ja, const int two_jb, const int two_jc);
static int lndelta_e (const int two_j1, const int two_j2, const int two_j3,
		      gsl_sf_result * lndelta);

/* This macro returns 1 if n is even, -1 if n is odd. */
#define PHASE(n) (GSL_IS_ODD(n) ? -1.0 : 1.0)

static int
istriangle (const int two_ja, const int two_jb, const int two_jc)
     /* Returns 0 if the triangle condition is met, and the sum of twice the
        angular momenta is even. Arguments are twice the value of the angular
        momenta. */
{
  if ((two_jc <= two_ja + two_jb) && (two_jc >= abs (two_ja - two_jb))
      && (GSL_IS_EVEN (two_ja + two_jb + two_jc)))
    return 1;
  else
    return 0;
}

static int
lndelta_e (const int two_j1, const int two_j2,
           const int two_j3, gsl_sf_result * lndelta)
     /* Calculates the natural log of the delta function - Zare Eq. A-2. Note:
        values returned from gsl_sf_ln_fact_e are always positive. */
{
  gsl_sf_result a1, a2, a3, a4;
  int status;

  status = gsl_sf_lnfact_e ((two_j1 + two_j2 - two_j3) / 2, &a1);
  status += gsl_sf_lnfact_e ((two_j1 - two_j2 + two_j3) / 2, &a2);
  status += gsl_sf_lnfact_e ((-two_j1 + two_j2 + two_j3) / 2, &a3);
  status += gsl_sf_lnfact_e ((two_j1 + two_j2 + two_j3) / 2 + 1, &a4);

  if (status != GSL_SUCCESS)
    OVERFLOW_ERROR (lndelta);

  lndelta->val = 0.5 * (a1.val + a2.val + a3.val - a4.val);
  lndelta->err = 0.5 * (a1.err + a2.err + a3.err + a4.err);
  lndelta->err += 2.0 * GSL_DBL_EPSILON * (a1.val + a2.val + a3.val + a4.val);

  return GSL_SUCCESS;
}

int
gsl_sf_wigner_3j_e (const int two_j1, const int two_j2, const int two_j3,
                    const int two_m1, const int two_m2, const int two_m3,
                    gsl_sf_result * result)
     /* Returns the value of the wigner 3j symbol - Zare Eq. A-1. Note that in
        Zare's book, this equation contains typos. */
{
  int v, vmin, vmax, t1, t2, t3, t4, t5, status;
  gsl_sf_result a, a1, a2, a3, a4, a5, a6, a7, a8;

  if ((two_j1 < 0) || (two_j2 < 0) || (two_j3 < 0) ||
      (abs (two_m1) > two_j1) || (abs (two_m2) > two_j2) ||
      (abs (two_m3) > two_j3) || GSL_IS_ODD (two_j1 + two_m1) ||
      GSL_IS_ODD (two_j2 + two_m2) || GSL_IS_ODD (two_j3 + two_m3))
    DOMAIN_ERROR (result);

  result->val = 0.0;
  result->err = 0.0;

  if ((!istriangle (two_j1, two_j2, two_j3)) ||
      (two_m1 + two_m2 + two_m3) != 0)
    return GSL_SUCCESS;

  t1 = (two_j1 + two_j2 - two_j3) / 2;
  t2 = (two_j1 - two_m1) / 2;
  t3 = (two_j2 + two_m2) / 2;
  t4 = (two_j3 - two_j2 + two_m1) / 2;
  t5 = (two_j3 - two_j1 - two_m2) / 2;

  vmin = (t4 < t5) ? -t4 : -t5;
  vmin = (vmin < 0) ? 0 : vmin;

  vmax = t1;
  vmax = (t2 < vmax) ? t2 : vmax;
  vmax = (t3 < vmax) ? t3 : vmax;

  if (vmin > vmax)
    return GSL_SUCCESS;

  status = gsl_sf_lnfact_e ((two_j1 + two_m1) / 2, &a1);
  status += gsl_sf_lnfact_e (t2, &a2);
  status += gsl_sf_lnfact_e (t3, &a3);
  status += gsl_sf_lnfact_e ((two_j2 - two_m2) / 2, &a4);
  status += gsl_sf_lnfact_e ((two_j3 + two_m3) / 2, &a5);
  status += gsl_sf_lnfact_e ((two_j3 - two_m3) / 2, &a6);
  status += lndelta_e (two_j1, two_j2, two_j3, &a7);

  if (status != GSL_SUCCESS)
    OVERFLOW_ERROR (result);

  a.val =
    0.5 * (a1.val + a2.val + a3.val + a4.val + a5.val + a6.val) + a7.val;
  a.err =
    0.5 * (a1.err + a2.err + a3.err + a4.err + a5.err + a6.err) + a7.err;
  a.err += 2.0 * GSL_DBL_EPSILON * a.val;

  for (v = vmin; v <= vmax; v++)
    {
      status += gsl_sf_lnfact_e (v, &a1);
      status += gsl_sf_lnfact_e (t1 - v, &a2);
      status += gsl_sf_lnfact_e (t2 - v, &a3);
      status += gsl_sf_lnfact_e (t3 - v, &a4);
      status += gsl_sf_lnfact_e (t4 + v, &a5);
      status += gsl_sf_lnfact_e (t5 + v, &a6);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      a7.val = a.val - (a1.val + a2.val + a3.val + a4.val + a5.val + a6.val);
      a7.err = (a.err + a1.err + a2.err + a3.err + a4.err + a5.err + a6.err);
      a7.err += 2.0 * GSL_DBL_EPSILON *
        (a.val + a1.val + a2.val + a3.val + a4.val + a5.val + a6.val);

      status += gsl_sf_exp_err_e (a7.val, a7.err, &a8);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      result->val += PHASE (v) * a8.val;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs (result->val);
    }

  result->val *= PHASE ((two_j1 - two_j2 - two_m3) / 2);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs (result->val);

  return GSL_SUCCESS;
}

int
gsl_sf_wigner_6j_e (const int two_j1, const int two_j2, const int two_j3,
                    const int two_j4, const int two_j5, const int two_j6,
                    gsl_sf_result * result)
     /* Returns the value of the 6j symbol - Zare Eq. A-4 */
{
  int v, vmin, vmax, t1, t2, t3, t4, t5, t6, t7, status;
  gsl_sf_result a, a1, a2, a3, a4;

  if ((two_j1 < 0) || (two_j2 < 0) || (two_j3 < 0) ||
      (two_j4 < 0) || (two_j5 < 0) || (two_j6 < 0))
    DOMAIN_ERROR (result);

  result->val = 0.0;
  result->err = 0.0;

  if (!istriangle (two_j1, two_j2, two_j3) ||
      !istriangle (two_j1, two_j5, two_j6) ||
      !istriangle (two_j4, two_j2, two_j6) ||
      !istriangle (two_j4, two_j5, two_j3))
    return GSL_SUCCESS;

  t1 = (two_j1 + two_j2 + two_j3) / 2;
  t2 = (two_j1 + two_j5 + two_j6) / 2;
  t3 = (two_j4 + two_j2 + two_j6) / 2;
  t4 = (two_j4 + two_j5 + two_j3) / 2;
  t5 = (two_j1 + two_j2 + two_j4 + two_j5) / 2;
  t6 = (two_j2 + two_j3 + two_j5 + two_j6) / 2;
  t7 = (two_j3 + two_j1 + two_j6 + two_j4) / 2;

  vmin = t1;
  vmin = (t2 > vmin) ? t2 : vmin;
  vmin = (t3 > vmin) ? t3 : vmin;
  vmin = (t4 > vmin) ? t4 : vmin;

  vmax = t5;
  vmax = (t6 < vmax) ? t6 : vmax;
  vmax = (t7 < vmax) ? t7 : vmax;

  if (vmin > vmax)
    return GSL_SUCCESS;

  status = lndelta_e (two_j1, two_j2, two_j3, &a1);
  status += lndelta_e (two_j1, two_j5, two_j6, &a2);
  status += lndelta_e (two_j4, two_j2, two_j6, &a3);
  status += lndelta_e (two_j4, two_j5, two_j3, &a4);

  if (status != GSL_SUCCESS)
    OVERFLOW_ERROR (result);

  a.val = a1.val + a2.val + a3.val + a4.val;
  a.err = a1.err + a2.err + a3.err + a4.err;
  a.err += 2.0 * GSL_DBL_EPSILON * a.val;

  for (v = vmin; v <= vmax; v++)
    {
      gsl_sf_result b1, b2, b3, b4, b5, b6, b7, b8, b9, b10;
      status += gsl_sf_lnfact_e (v + 1, &b1);
      status += gsl_sf_lnfact_e (v - t1, &b2);
      status += gsl_sf_lnfact_e (v - t2, &b3);
      status += gsl_sf_lnfact_e (v - t3, &b4);
      status += gsl_sf_lnfact_e (v - t4, &b5);
      status += gsl_sf_lnfact_e (t5 - v, &b6);
      status += gsl_sf_lnfact_e (t6 - v, &b7);
      status += gsl_sf_lnfact_e (t7 - v, &b8);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      b9.val = a.val + b1.val - b2.val - b3.val - b4.val - b5.val - b6.val -
        b7.val - b8.val;
      b9.err = a.err + b1.err + b2.err + b3.err + b4.err + b5.err + b6.err +
        b7.err + b8.err;
      b9.err += 2.0 * GSL_DBL_EPSILON * (a.val + b1.val + b2.val + b3.val +
                                         b4.val + b5.val + b6.val + b7.val +
                                         b8.val);

      status += gsl_sf_exp_err_e (b9.val, b9.err, &b10);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      result->val += b10.val * PHASE (v);
      result->err += b10.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs (result->val);
    }

  return GSL_SUCCESS;
}

int
gsl_sf_wigner_9j_e (const int two_j1, const int two_j2, const int two_j3,
                    const int two_j4, const int two_j5, const int two_j6,
                    const int two_j7, const int two_j8, const int two_j9,
                    gsl_sf_result * result)
    /* Returns the value of 9j symbol - Zare Eq. 4.24 */
{
  int two_k, two_kmin, two_kmax, status = GSL_SUCCESS;

  if ((two_j1 < 0) || (two_j2 < 0) || (two_j3 < 0) ||
      (two_j4 < 0) || (two_j5 < 0) || (two_j6 < 0) ||
      (two_j7 < 0) || (two_j8 < 0) || (two_j9 < 0))
    DOMAIN_ERROR (result);

  result->val = 0.0;
  result->err = 0.0;

  if (!istriangle (two_j1, two_j2, two_j3) ||
      !istriangle (two_j4, two_j5, two_j6) ||
      !istriangle (two_j7, two_j8, two_j9) ||
      !istriangle (two_j1, two_j4, two_j7) ||
      !istriangle (two_j2, two_j5, two_j8) ||
      !istriangle (two_j3, two_j6, two_j9))
    return GSL_SUCCESS;

  two_kmin = abs (two_j1 - two_j9);
  two_kmin =
    (abs (two_j8 - two_j4) < two_kmin) ? abs (two_j8 - two_j4) : two_kmin;
  two_kmin =
    (abs (two_j6 - two_j2) < two_kmin) ? abs (two_j6 - two_j2) : two_kmin;

  two_kmax = (two_j1 + two_j9);
  two_kmax = (two_j8 + two_j4 > two_kmax) ? two_j8 + two_j4 : two_kmax;
  two_kmax = (two_j6 + two_j2 > two_kmax) ? two_j6 + two_j2 : two_kmax;

  if (two_kmax < two_kmin)
    {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }

  for (two_k = two_kmin; two_k <= two_kmax; two_k += 2)
    {
      gsl_sf_result a1, a2, a3;
      status += gsl_sf_wigner_6j_e (two_j1, two_j4, two_j7,
                                    two_j8, two_j9, two_k, &a1);
      status +=
        gsl_sf_wigner_6j_e (two_j2, two_j5, two_j8,
                            two_j4, two_k, two_j6, &a2);
      status +=
        gsl_sf_wigner_6j_e (two_j3, two_j6, two_j9,
                            two_k, two_j1, two_j2, &a3);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      result->val += (two_k + 1.0) * PHASE (two_k) * a1.val * a2.val * a3.val;

      result->err += result->val * (fabs (a1.err / a1.val) +
                                    fabs (a2.err / a2.val) +
                                    fabs (a3.err / a3.val));
      result->err += 2.0 * GSL_DBL_EPSILON * fabs (result->val);
    }
  return GSL_SUCCESS;
}

int
gsl_sf_wigner_drot_e (const int two_j, const int two_m1, const int two_m2,
                      const double theta, gsl_sf_result * result)
     /* Returns the value of the reduced rotation matrix - Zare Eq. 3.57. (which
        contains typos).  Angle theta is in radians. */
{
  int v, vmin, vmax, t1, t2, t3, t4, tv, status = GSL_SUCCESS;
  gsl_sf_result a, a1, a2, a3, a4, st, ct;

  if ((two_j) < 0 || abs (two_m1) > two_j || abs (two_m2) > two_j ||
      (GSL_IS_ODD (two_j + two_m1)) || (GSL_IS_ODD (two_j + two_m2)))
    DOMAIN_ERROR (result);

  result->val = 0.0;
  result->err = 0.0;

  t1 = (two_j + two_m1) / 2;
  t2 = (two_j - two_m2) / 2;
  t3 = (two_m2 - two_m1) / 2;
  t4 = (2 * two_j + two_m1 - two_m2) / 2;

  vmin = (-t3 < 0) ? 0 : -t3;
  vmax = (t2 < t1) ? t2 : t1;

  if (vmin > vmax)
    return GSL_SUCCESS;

  status += gsl_sf_cos_e (theta * 0.5, &ct);
  status += gsl_sf_sin_e (theta * 0.5, &st);
  if (status != GSL_SUCCESS)
    OVERFLOW_ERROR (result);

  status += gsl_sf_lnfact_e (t1, &a1);
  status += gsl_sf_lnfact_e (t2, &a2);
  status += gsl_sf_lnfact_e ((two_j - two_m1) / 2, &a3);
  status += gsl_sf_lnfact_e ((two_j + two_m2) / 2, &a4);
  if (status != GSL_SUCCESS)
    OVERFLOW_ERROR (result);

  a.val = 0.5 * (a1.val + a2.val + a3.val + a4.val);
  a.err = 0.5 * (a1.err + a2.err + a3.err + a4.err);
  a.err += 2.0 * GSL_DBL_EPSILON * a.val;

  for (v = vmin; v <= vmax; v++)
    {
      gsl_sf_result b1, b2, b3, b4, b5, b6, b7, b8;
      int i1, i2;

      status += gsl_sf_lnfact_e (t1 - v, &b1);
      status += gsl_sf_lnfact_e (t2 - v, &b2);
      status += gsl_sf_lnfact_e (t3 + v, &b3);
      status += gsl_sf_lnfact_e (v, &b4);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      b5.val = a.val - b1.val - b2.val - b3.val - b4.val;
      b5.err = a.err + b1.err + b2.err + b3.err + b4.err;
      b5.err +=
        2.0 * GSL_DBL_EPSILON * (a.val + b1.val + b2.val + b3.val + b4.val);

      status += gsl_sf_exp_err_e (b5.val, b5.err, &b6);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      tv = 2 * v;
      i1 = t4 - tv;
      i2 = t3 + tv;

      status += gsl_sf_pow_int_e (ct.val, i1, &b7);
      status += gsl_sf_pow_int_e (st.val, i2, &b8);
      if (status != GSL_SUCCESS)
        OVERFLOW_ERROR (result);

      /* Bolt on error in pow_int for the error in the input values */
      b7.err += i1 * fabs (ct.err);
      b8.err += i2 * fabs (st.err);

      result->val += PHASE (v) * b6.val * b7.val * b8.val;

      result->err += fabs (result->val) *
        (fabs (b6.err / b6.val) + fabs (b7.err / b7.val) +
         fabs (b8.err / b8.val));
      result->err += 2.0 * GSL_DBL_EPSILON * fabs (result->val);
    }

  return GSL_SUCCESS;
}

#include "eval.h"

double
gsl_sf_wigner_3j (const int two_j1, const int two_j2, const int two_j3,
                  const int two_m1, const int two_m2, const int two_m3)
{
  EVAL_RESULT (gsl_sf_wigner_3j_e (two_j1, two_j2, two_j3,
                                   two_m1, two_m2, two_m3, &result));
}

double
gsl_sf_wigner_6j (const int two_j1, const int two_j2, const int two_j3,
                  const int two_j4, const int two_j5, const int two_j6)
{
  EVAL_RESULT (gsl_sf_wigner_6j_e (two_j1, two_j2, two_j3,
                                   two_j4, two_j5, two_j6, &result));
}

double
gsl_sf_wigner_9j (const int two_j1, const int two_j2, const int two_j3,
                  const int two_j4, const int two_j5, const int two_j6,
                  const int two_j7, const int two_j8, const int two_j9)
{
  EVAL_RESULT (gsl_sf_wigner_9j_e (two_j1, two_j2, two_j3,
                                   two_j4, two_j5, two_j6,
                                   two_j7, two_j8, two_j9, &result));
}

double
gsl_sf_wigner_drot (const int two_j, const int two_m1, const int two_m2,
                    const double theta)
{
  EVAL_RESULT (gsl_sf_wigner_drot_e (two_j, two_m1, two_m2, theta, &result));
}

#undef PHASE

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* matrix/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#define M 53
#define N 107

int status = 0;

#ifndef DESC
#define DESC ""
#endif

#define BASE_GSL_COMPLEX_LONG
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_LONG

#define BASE_GSL_COMPLEX
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX

#define BASE_GSL_COMPLEX_FLOAT
#include "templates_on.h"
#include "test_complex_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_FLOAT

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "test_source.c"
#include "templates_off.h"
#undef  BASE_CHAR

void my_error_handler (const char *reason, const char *file,
		       int line, int err);

int
main (void)
{
  gsl_ieee_env_setup ();

  test_func ();
  test_float_func ();
  test_long_double_func ();
  test_ulong_func ();
  test_long_func ();
  test_uint_func ();
  test_int_func ();
  test_ushort_func ();
  test_short_func ();
  test_uchar_func ();
  test_char_func ();
  test_complex_func ();
  test_complex_float_func ();
  test_complex_long_double_func ();

  test_text ();
  test_float_text ();
#ifdef HAVE_PRINTF_LONGDOUBLE
  test_long_double_text ();
#endif
  test_ulong_text ();
  test_long_text ();
  test_uint_text ();
  test_int_text ();
  test_ushort_text ();
  test_short_text ();
  test_uchar_text ();
  test_char_text ();
  test_complex_text ();
  test_complex_float_text ();
#ifdef HAVE_PRINTF_LONGDOUBLE
  test_complex_long_double_text ();
#endif

  test_binary ();
  test_float_binary ();
  test_long_double_binary ();
  test_ulong_binary ();
  test_long_binary ();
  test_uint_binary ();
  test_int_binary ();
  test_ushort_binary ();
  test_short_binary ();
  test_uchar_binary ();
  test_char_binary ();
  test_complex_binary ();
  test_complex_float_binary ();
  test_complex_long_double_binary ();

  gsl_warnings_off = 1;

  gsl_set_error_handler (&my_error_handler);

  test_trap ();
  test_float_trap ();
  test_long_double_trap ();
  test_ulong_trap ();
  test_long_trap ();
  test_uint_trap ();
  test_int_trap ();
  test_ushort_trap ();
  test_short_trap ();
  test_uchar_trap ();
  test_char_trap ();
  test_complex_trap ();
  test_complex_float_trap ();
  test_complex_long_double_trap ();

  exit (gsl_test_summary ());
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
  status = 1;
}


text:
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>

#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_halfcomplex_float.h>

#define BASE_DOUBLE
#include "templates_on.h"
#include "hc_pass.h"
#include "hc_init.c"
#include "hc_main.c"
#include "hc_pass_2.c"
#include "hc_pass_3.c"
#include "hc_pass_4.c"
#include "hc_pass_5.c"
#include "hc_pass_n.c"
#include "hc_radix2.c"
#include "hc_unpack.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "hc_pass.h"
#include "hc_init.c"
#include "hc_main.c"
#include "hc_pass_2.c"
#include "hc_pass_3.c"
#include "hc_pass_4.c"
#include "hc_pass_5.c"
#include "hc_pass_n.c"
#include "hc_radix2.c"
#include "hc_unpack.c"
#include "templates_off.h"
#undef  BASE_FLOAT


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_neurotrace.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_stack_io.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"
#include "tz_stack_threshold.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];

  sprintf(file_path, "../data/%s.tif", neuron_name);
  Stack *signal = Read_Stack(file_path);
  Translate_Stack(signal, GREY, 1);
  Stack *stack = Translate_Stack(signal, COLOR, 0);

  sprintf(file_path, "../data/%s/soma.tif", neuron_name);
  Stack *label = Read_Stack(file_path);
  Stack *label2 = Copy_Stack(label);

  Stack_Threshold(label, 2);
  Stack_Binarize(label);

  Stack_Threshold(label2, 1);
  Stack_Binarize(label2);
  Stack_Xor(label, label2, label2);

  Print_Stack_Info(signal);

  Stack_Label_Color(stack, label, 0.0, 1.0, signal);
  Stack_Label_Color(stack, label2, 4.0, 1.0, signal);

  int i;
  for (i = 0; i < 3; i++) {
    Stack_Channel_Extraction(stack, i, signal);
    sprintf(file_path, "../data/%s/label_%d.tif", neuron_name, i);
    Write_Stack(file_path, signal);
  }

  Kill_Stack(signal);
  Kill_Stack(stack);
  Kill_Stack(label);
  Kill_Stack(label2);

  sprintf(file_path, "../data/%s/traced.tif", neuron_name);
  Stack *canvas = Read_Stack(file_path);
  sprintf(file_path, "../data/%s/soma.tif", neuron_name);
  Stack *soma = Read_Stack(file_path);
  Stack_Blend_Mc_L(canvas, soma, 2, 4.0);
  Stack_Blend_Mc_L(canvas, soma, 3, 1.0);
  sprintf(file_path, "../data/%s/label_color.tif", neuron_name);
  Write_Stack(file_path, canvas);

  return 0;
}


text:
/* blas/gsl_blas.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
 * Author:  G. Jungman
 */
#ifndef __GSL_BLAS_H__
#define __GSL_BLAS_H__

#include <gsl/vector/gsl_vector.h>
#include <gsl/matrix/gsl_matrix.h>

#include <gsl/blas/gsl_blas_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/**
	MODIFIED BY STELIAN:
		- only keep a subset of the methods that I think are important
		- link to an external blas library
*/

/* ========================================================================
 * Level 1
 * ========================================================================
 */

/**
*     forms the dot product of two vectors.
*/
GSL_FUNC int gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double * result);

/**
*  DNRM2 returns the euclidean norm of a vector via the function
*  name, so that
*
*     DNRM2 := sqrt( x'*x )
*/
GSL_FUNC double gsl_blas_dnrm2  (const gsl_vector * X);

/**
*     interchanges two vectors.
*/
GSL_FUNC int  gsl_blas_dswap (gsl_vector * X, gsl_vector * Y);

/**
*     copies a vector, x, to a vector, y.
*/
GSL_FUNC int  gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y);

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

/**
*  DGEMV  performs one of the matrix-vector operations
*
*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
*
*  where alpha and beta are scalars, x and y are vectors and A is an
*  m by n matrix.
*/
GSL_FUNC int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A, const gsl_vector * X, double beta, gsl_vector * Y);

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/**
*  DGEMM  performs one of the matrix-matrix operations
*
*     C := alpha*op( A )*op( B ) + beta*C,
*
*  where  op( X ) is one of
*
*     op( X ) = X   or   op( X ) = X',
*
*  alpha and beta are scalars, and A, B and C are matrices, with op( A )
*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
*/
GSL_FUNC int  gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, double alpha, const gsl_matrix * A, const gsl_matrix * B, double beta, gsl_matrix * C);

__END_DECLS

#endif /* __GSL_BLAS_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SIM_MATH_ODE_SOLVER_H_
#define SIM_MATH_ODE_SOLVER_H_

#include <stdint.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

#include <vector>

namespace sim {

// Enumeration for return codes.
enum class OdeSolverStatus { kSuccess, kError, kOutOfPhysics };

// Interface representing a system of ordinary differential equations.
class OdeSystem {
 public:
  virtual ~OdeSystem() {}

  // Get the number of states for the ODE.
  virtual int32_t num_states() const = 0;

  // Calculate the state derivative for a given time and state.
  virtual void CalcDerivatives(double t, const std::vector<double> &x,
                               std::vector<double> *dx) const = 0;
};

// Interface for ODE solvers.
class OdeSolver {
 public:
  virtual ~OdeSolver() {}

  // Integrate the ODE with initial condition x0 at time t0 until time tf.
  //
  // Args:
  //   t0: Initial time.
  //   tf: Final time (>= t0).
  //   x0: Initial state.
  //   t_int: Null pointer or output time at which the integrator terminated.
  //   x: Terminal state.  It should be safe to reuse x0 as x.
  //
  // Returns:
  //   kSuccess on success (indicating *t_int = tf) and kError otherwise.
  virtual OdeSolverStatus Integrate(double t0, double tf,
                                    const std::vector<double> &x0,
                                    double *t_int, std::vector<double> *x) = 0;
};

}  // namespace sim

#endif  // SIM_MATH_ODE_SOLVER_H_


text:
#ifndef __CCL_CONSTANTS_H_INCLUDED__
#define __CCL_CONSTANTS_H_INCLUDED__

#include <gsl/gsl_const_mksa.h>

//Spline types
#define A_SPLINE_TYPE gsl_interp_akima
#define K_SPLINE_TYPE gsl_interp_akima
#define L_SPLINE_TYPE gsl_interp_akima
#define M_SPLINE_TYPE gsl_interp_akima
#define D_SPLINE_TYPE gsl_interp_akima
#define PNL_SPLINE_TYPE gsl_interp2d_bicubic
#define PLIN_SPLINE_TYPE gsl_interp2d_bicubic
#define CORR_SPLINE_TYPE gsl_interp_akima

/** @file */

#ifndef M_PI
/**
 *  PI (in case it's not defined from math.h)
*/
#define M_PI 3.14159265358979323846
#endif

/**
 *  k pivot. These are in units of Mpc (no factor of h)
*/
#define K_PIVOT 0.05

/**
 * Lightspeed / H0 in units of Mpc/h (from CODATA 2014)
 */
#define CLIGHT_HMPC 2997.92458 //H0^-1 in Mpc/h

/**
 * Newton's gravitational constant in units of m^3/Kg/s^2
 */
//#define GNEWT 6.6738e-11  /(from PDG 2013) in m^3/Kg/s^2
//#define GNEWT 6.67428e-11 // CLASS VALUE
#define GNEWT 6.67408e-11 // from CODATA 2014

/**
 * Solar mass in units of kg (from GSL)
 */
//#define SOLAR_MASS GSL_CONST_MKSA_SOLAR_MASS
//#define SOLAR_MASS 1.9885e30 //(from PDG 2015) in Kg
#define SOLAR_MASS 1.9884754153381438e+30 //from IAU 2015
/**
 * Mpc to meters (from PDG 2013)
 */
#define MPC_TO_METER 3.08567758149e22

/**
 * pc to meters (from PDG 2013)
 */
#define PC_TO_METER 3.08567758149e16

/**
 * Rho critical in units of M_sun/h / (Mpc/h)^3
 */
#define RHO_CRITICAL ((3*100*100)/(8*M_PI*GNEWT)) * (1000*1000*MPC_TO_METER/SOLAR_MASS)

/**
 * Boltzmann constant in units of J/K
*/
//#define KBOLTZ  GSL_CONST_MKSA_BOLTZMANN
#define KBOLTZ 1.38064852e-23 //from CODATA 2014

/**
 * Stefan-Boltzmann constant in units of kg/s^3 / K^4
 */
//#define STBOLTZ GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT
#define STBOLTZ 5.670367e-8 //from CODATA 2014
/**
 * Planck's constant in units kg m^2 / s
 */
//#define HPLANCK  GSL_CONST_MKSA_PLANCKS_CONSTANT_H
#define HPLANCK 6.626070040e-34 //from CODATA 2014

/**
 * The speed of light in m/s
 */
//#define CLIGHT   GSL_CONST_MKSA_SPEED_OF_LIGHT
#define CLIGHT 299792458.0 //from CODATA 2014

/**
 * Electron volt to Joules convestion
 */
//#define EV_IN_J  GSL_CONST_MKSA_ELECTRON_VOLT
#define EV_IN_J 1.6021766208e-19  //from CODATA 2014

/**
 * Temperature of the CMB in K
 */
#define TCMB 2.725
//#define TCMB 2.7255 // CLASS value

/**
 * T_ncdm, as taken from CLASS, explanatory.ini
 */
#define TNCDM 0.71611

/**
 * neutrino mass splitting differences
 * See Lesgourgues and Pastor, 2012 for these values.
 * Adv. High Energy Phys. 2012 (2012) 608515,
 * arXiv:1212.6154, page 13
*/
#define DELTAM12_sq 7.62E-5
#define DELTAM13_sq_pos 2.55E-3
#define DELTAM13_sq_neg -2.43E-3

//Precision parameters
/**
 * Default relative precision if not otherwise specified
 */
#define GSL_EPSREL 1E-4

/**
 * Default number of iterations for integration and root-finding if not otherwise
 * specified
 */
#define GSL_N_ITERATION 1000

/**
 * Default number of Gauss-Kronrod points in QAG integration if not otherwise
 * specified
 */
#define GSL_INTEGRATION_GAUSS_KRONROD_POINTS GSL_INTEG_GAUSS41

/**
 * Absolute precision in neutrino root finding
 */
#define GSL_EPSABS_NU 1E-7

/**
 * Relative precision in neutrino root finding
 */
#define GSL_EPSREL_NU 1E-7

/**
 * Number of iterations for neutrino root finding
 */
#define GSL_N_ITERATION_NU 1000

/**
 * Relative precision in sigma_R calculations
 */
#define GSL_EPSREL_SIGMAR 1E-5

/**
 * Relative precision in distance calculations
 */
#define GSL_EPSREL_DIST 1E-6

/**
 * Relative precision in growth calculations
 */
#define GSL_EPSREL_GROWTH 1E-6

/**
 * Relative precision in dNdz calculations
 */
#define GSL_EPSREL_DNDZ 1E-6

/**
 * Absolute precision in growth calculations
 */
#define EPS_SCALEFAC_GROWTH 1E-6

#endif

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
/* movstat/madacc.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.
 */

/*
 * This module contains routines for computing the median absolute deviation (MAD)
 * of a moving fixed-sized window.
 */

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

typedef double madacc_type_t;
typedef madacc_type_t ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  size_t n;                        /* window size */
  const gsl_movstat_accum *medacc; /* median accumulator */
  void *median_state;              /* median accumulator workspace */
  ringbuf *rbuf;                   /* ring buffer storing current window, size n */
  madacc_type_t *work;             /* workspace, size n */
} madacc_state_t;

static size_t madacc_size(const size_t n);
static int madacc_init(const size_t n, void * vstate);
static int madacc_insert(const madacc_type_t x, void * vstate);
static int madacc_delete(void * vstate);
static int madacc_medmad(void * params, madacc_type_t * result, const void * vstate);

static size_t
madacc_size(const size_t n)
{
  size_t size = 0;
  const gsl_movstat_accum * acc = gsl_movstat_accum_median;

  size += sizeof(madacc_state_t);
  size += (acc->size)(n);            /* median accumulator */
  size += ringbuf_size(n);           /* rbuf */
  size += n * sizeof(madacc_type_t); /* work */

  return size;
}

static int
madacc_init(const size_t n, void * vstate)
{
  madacc_state_t * state = (madacc_state_t *) vstate;

  state->n = n;
  state->medacc = gsl_movstat_accum_median;

  state->median_state = (void *) ((unsigned char *) vstate + sizeof(madacc_state_t));
  state->rbuf = (ringbuf *) ((unsigned char *) state->median_state + (state->medacc->size)(n));
  state->work = (madacc_type_t *) ((unsigned char *) state->rbuf + ringbuf_size(n));

  (state->medacc->init)(n, state->median_state);
  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
madacc_insert(const madacc_type_t x, void * vstate)
{
  madacc_state_t * state = (madacc_state_t *) vstate;

  /* insert element into median accumulator */
  (state->medacc->insert)(x, state->median_state);

  /* insert element into ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
madacc_delete(void * vstate)
{
  madacc_state_t * state = (madacc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    ringbuf_pop_back(state->rbuf);

  return GSL_SUCCESS;
}

static int
madacc_medmad(void * params, madacc_type_t * result, const void * vstate)
{
  const madacc_state_t * state = (const madacc_state_t *) vstate;

  if (ringbuf_is_empty(state->rbuf))
    {
      GSL_ERROR("no samples yet added to workspace", GSL_EINVAL);
    }
  else
    {
      int status;
      const double scale = *(double *) params;
      const int n = ringbuf_n(state->rbuf);
      int i;
      double median, mad;

      /* compute median of current window */
      status = (state->medacc->get)(NULL, &median, state->median_state);
      if (status)
        return status;

      /* compute current window absolute deviations from median */
      for (i = 0; i < n; ++i)
        {
          double xi = state->rbuf->array[(state->rbuf->head + i) % state->rbuf->size];
          state->work[i] = fabs(xi - median);
        }

      /* compute MAD of current window */
      mad = gsl_stats_median(state->work, 1, n);

      result[0] = median;
      result[1] = scale * mad;

      return GSL_SUCCESS;
    }
}

static const gsl_movstat_accum mad_accum_type =
{
  madacc_size,
  madacc_init,
  madacc_insert,
  NULL, /*XXX*/
  madacc_medmad
};

const gsl_movstat_accum *gsl_movstat_accum_mad = &mad_accum_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

/*
 * Spline the linear power spectrum for mu-Sigma MG cosmologies.
 * @param cosmo Cosmological parameters
 ^ @param psp The linear power spectrum to spline.
 * @param status, integer indicating the status
 */
void ccl_cosmology_spline_linpower_musigma(ccl_cosmology* cosmo, ccl_f2d_t *psp, int rescaled_mg_flag, int* status) {
  double kmin, kmax, ndecades, amin, amax, ic, sigma8, log_sigma8;
  int nk, na, s;
  double *lk = NULL, *aa = NULL, *lpk_ln = NULL, *lpk_nl = NULL;
  double norm_pk;
  double *mnu_list = NULL;
  ccl_parameters params_GR;
  params_GR.m_nu = NULL;
  params_GR.z_mgrowth = NULL;
  params_GR.df_mgrowth = NULL;
  ccl_cosmology * cosmo_GR = NULL;
  double *D_mu = NULL;
  double *D_GR = NULL;

  if (*status == 0) {
    //calculations done - now allocate CCL splines
    kmin = 2*exp(psp->lkmin);
    kmax = cosmo->spline_params.K_MAX_SPLINE;
    //Compute nk from number of decades and N_K = # k per decade
    ndecades = log10(kmax) - log10(kmin);
    nk = (int)ceil(ndecades*cosmo->spline_params.N_K);
    amin = cosmo->spline_params.A_SPLINE_MINLOG_PK;
    amax = cosmo->spline_params.A_SPLINE_MAX;
    na = cosmo->spline_params.A_SPLINE_NA_PK+cosmo->spline_params.A_SPLINE_NLOG_PK-1;

    // The lk array is initially k, but will later
    // be overwritten with log(k)
    lk = ccl_log_spacing(kmin, kmax, nk);
    if (lk == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
          cosmo,
          "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
    }
  }

  if (*status == 0) {
    aa = ccl_linlog_spacing(
        amin, cosmo->spline_params.A_SPLINE_MIN_PK,
        amax, cosmo->spline_params.A_SPLINE_NLOG_PK,
        cosmo->spline_params.A_SPLINE_NA_PK);
    if (aa == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
    }
  }

  if (*status == 0) {
    lpk_ln = malloc(nk * na * sizeof(double));
    if (lpk_ln == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
    }
  }

  if (*status == 0) {
    lpk_nl = malloc(nk * na * sizeof(double));
    if(lpk_nl == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
    }
  }

  if (*status == 0) {
    // After this loop lk will contain log(k),
    // lpk_ln will contain log(P_lin), all in Mpc, not Mpc/h units!
    double psout_l;
    s = 0;

    // If scale-independent mu / Sigma modified gravity is in use
    // and mu ! = 0 : get the unnormalized growth factor in MG and for
    // corresponding GR case, to rescale CLASS power spectrum
    if (fabs(cosmo->params.mu_0) > 1e-14) {
      // Set up another cosmology which is exactly the same as the
      // current one but with mu_0 and Sigma_0=0, for scaling P(k)

      // Get a list of the three neutrino masses already calculated
      mnu_list = malloc(3*sizeof(double));
      if (mnu_list == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
      }

      if (*status == 0) {
        for (int i=0; i< cosmo->params.N_nu_mass; i=i+1) {
          mnu_list[i] = cosmo->params.m_nu[i];
        }
        if (cosmo->params.N_nu_mass < 3) {
          for (int j=cosmo->params.N_nu_mass; j<3; j=j+1) {
            mnu_list[j] = 0.;
          }
        }

        if (isfinite(cosmo->params.A_s)) {
          norm_pk = cosmo->params.A_s;
        }
        else if (isfinite(cosmo->params.sigma8)) {
          norm_pk = cosmo->params.sigma8;
        }
        else {
          *status = CCL_ERROR_PARAMETERS;
          strcpy(
            cosmo->status_message,
            "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): neither A_s nor sigma8 defined.\n");
        }
      }

      if (*status == 0) {
        params_GR = ccl_parameters_create(
          cosmo->params.Omega_c, cosmo->params.Omega_b, cosmo->params.Omega_k,
          cosmo->params.Neff, mnu_list, cosmo->params.N_nu_mass,
          cosmo->params.w0, cosmo->params.wa, cosmo->params.h,
          norm_pk, cosmo->params.n_s,
          cosmo->params.bcm_log10Mc, cosmo->params.bcm_etab,
          cosmo->params.bcm_ks, 0., 0., cosmo->params.nz_mgrowth,
          cosmo->params.z_mgrowth, cosmo->params.df_mgrowth, status);

        if (*status) {
          *status = CCL_ERROR_PARAMETERS;
          strcpy(
            cosmo->status_message,
            "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): could not make MG params.\n");
        }
      }

      if (*status == 0) {
        cosmo_GR = ccl_cosmology_create(params_GR, cosmo->config);
        D_mu = malloc(na * sizeof(double));
        D_GR = malloc(na * sizeof(double));

        if (cosmo_GR == NULL || D_mu == NULL || D_GR == NULL) {
          *status = CCL_ERROR_MEMORY;
          ccl_cosmology_set_status_message(
            cosmo,
            "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): memory allocation\n");
        }
      }

      if (*status == 0) {
        ccl_cosmology_compute_growth(cosmo_GR, status);

        if (*status) {
          *status = CCL_ERROR_PARAMETERS;
          strcpy(
            cosmo->status_message,
            "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): could not init GR growth.\n");
        }
      }

      if (*status == 0) {
        for (int i=0; i<na; i++) {
          D_mu[i] = ccl_growth_factor_unnorm(cosmo, aa[i], status);
          D_GR[i] = ccl_growth_factor_unnorm(cosmo_GR, aa[i], status);
        }

        if (*status) {
          *status = CCL_ERROR_PARAMETERS;
          strcpy(
            cosmo->status_message,
            "ccl_power.c: ccl_cosmology_spline_linpower_musigma(): could not make MG and GR growth.\n");
        }
      }

      if (*status == 0) {
        for (int i=0; i<nk; i++) {
          lk[i] = log(lk[i]);
          for (int j = 0; j < na; j++) {
            //The 2D interpolation routines access the function values pk_{k_ia_j} with the following ordering:
            //pk_ij = pk[j*N_k + i]
            //with i = 0,...,N_k-1 and j = 0,...,N_a-1.
            psout_l = ccl_f2d_t_eval(psp, lk[i], aa[j], cosmo, status);
	    if (rescaled_mg_flag == 0) {
            lpk_ln[j*nk+i] = log(psout_l) ;
	    }
	    else {
            lpk_ln[j*nk+i] = log(psout_l) + 2 * log(D_mu[j]) - 2 * log(D_GR[j]);
	    }

          }
        }
      }
    }
    else {
      // This is the normal GR case.
      for (int i=0; i<nk; i++) {
        lk[i] = log(lk[i]);
        for (int j = 0; j<na; j++) {
          //The 2D interpolation routines access the function values pk_{k_ia_j} with the following ordering:
          //pk_ij = pk[j*N_k + i]
          //with i = 0,...,N_k-1 and j = 0,...,N_a-1.
          psout_l = ccl_f2d_t_eval(psp, lk[i], aa[j], cosmo, status);
          lpk_ln[j*nk+i] = log(psout_l);
        }
      }
    }
  }

  if (*status == 0) {
    cosmo->data.p_lin = ccl_f2d_t_new(
      na, aa, nk, lk, lpk_ln, NULL, NULL, 0,
      1, 2, ccl_f2d_cclgrowth, 1, NULL, 0, 2,
      ccl_f2d_3,status);
  }

  // if desried, renomalize to a given sigma8
  if (isfinite(cosmo->params.sigma8) && (!isfinite(cosmo->params.A_s))) {
    if (*status == 0) {
      cosmo->computed_linear_power = true;
      sigma8 = ccl_sigma8(cosmo, status);
      cosmo->computed_linear_power = false;
    }

    if (*status == 0) {
      // Calculate normalization factor using computed value of sigma8, then
      // recompute P(k, a) using this normalization
      log_sigma8 = 2*(log(cosmo->params.sigma8) - log(sigma8));
      for(int j = 0; j<na*nk; j++)
        lpk_ln[j] += log_sigma8;
    }

    if (*status == 0) {
      // Free the previous P(k,a) spline, and allocate a new one to store the
      // properly-normalized P(k,a)
      ccl_f2d_t_free(cosmo->data.p_lin);
      cosmo->data.p_lin = ccl_f2d_t_new(
        na, aa, nk, lk, lpk_ln, NULL, NULL, 0,
        1, 2, ccl_f2d_cclgrowth, 1, NULL, 0, 2,
        ccl_f2d_3,status);
    }
  }

  free(D_mu);
  free(D_GR);
  free(mnu_list);
  ccl_parameters_free(&params_GR);
  ccl_cosmology_free(cosmo_GR);
  free(lk);
  free(aa);
  free(lpk_nl);
  free(lpk_ln);
}

max_stars_repo_name: timothydmorton/CCL

max_issues_repo_name: timothydmorton/CCL

max_forks_repo_name: timothydmorton/CCL


text:
#include <string.h>

#include <mpi.h>
#include <pfft.h>
#include <gsl/gsl_rng.h>

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>
#include "pmpfft.h"

#define HAS(a, b) ((a & b) != 0)

static void
pack_any(FastPMStore * p, ptrdiff_t index, int ci, void * packed)
{
    size_t elsize = p->_column_info[ci].elsize;
    memcpy(packed, p->columns[ci] + index * elsize, elsize);
}

static void
unpack_any(FastPMStore * p, ptrdiff_t index, int ci, void * packed)
{
    size_t elsize = p->_column_info[ci].elsize;
    memcpy(p->columns[ci] + index * elsize, packed, elsize);
}

void
FastPMReduceOverwriteAny(FastPMStore * p, ptrdiff_t index, int ci, void * packed, void * userdata)
{
    size_t elsize = p->_column_info[ci].elsize;

    memcpy(p->columns[ci] + index * elsize, packed, elsize);
}

void
FastPMReduceAddFloat(FastPMStore * src, ptrdiff_t isrc, FastPMStore * dest, ptrdiff_t idest, int ci, void * userdata)
{
    size_t nmemb = dest->_column_info[ci].nmemb ;

    size_t elsize = dest->_column_info[ci].elsize;

    float * pdest = (float*) (dest->columns[ci] + idest * elsize);
    float * psrc  = (float*) (src ->columns[ci] + isrc * elsize);

    int d;
    for(d = 0; d < nmemb; d ++) {
        pdest[d] += psrc[d];
    }
}

static double
to_double_f4 (FastPMStore * p, ptrdiff_t index, int ci, int memb)
{
    size_t nmemb = p->_column_info[ci].nmemb ;
    if(memb > nmemb) {
        fastpm_raise(-1, "memb %d greater than nmemb %d", memb, nmemb);
    }
    size_t elsize = p->_column_info[ci].elsize;

    float * ptr = (float*) (p->columns[ci] + index * elsize);

    return ptr[memb];
}

static double
to_double_f8 (FastPMStore * p, ptrdiff_t index, int ci, int memb)
{
    size_t nmemb = p->_column_info[ci].nmemb ;
    if(memb > nmemb) {
        fastpm_raise(-1, "memb %d greater than nmemb %d", memb, nmemb);
    }
    size_t elsize = p->_column_info[ci].elsize;

    double * ptr = (double*) (p->columns[ci] + index * elsize);

    return ptr[memb];
}

static void
from_double_f4 (FastPMStore * p, ptrdiff_t index, int ci, int memb, const double value)
{
    size_t nmemb = p->_column_info[ci].nmemb ;
    if(memb > nmemb) {
        fastpm_raise(-1, "memb %d greater than nmemb %d", memb, nmemb);
    }
    size_t elsize = p->_column_info[ci].elsize;

    float * ptr = (float*) (p->columns[ci] + index * elsize);

    ptr[memb] = value;
}

const char *
fastpm_species_get_name(enum FastPMSpecies species)
{
    switch(species) {
        case FASTPM_SPECIES_BARYON:
            return "0";
        case FASTPM_SPECIES_CDM:
            return "1";
        case FASTPM_SPECIES_NCDM:
            return "2";
    }
    return "UNKNOWN";
}
void fastpm_store_get_position(FastPMStore * p, ptrdiff_t index, double pos[3])
{
    pos[0] = p->x[index][0];
    pos[1] = p->x[index][1];
    pos[2] = p->x[index][2];
}

void fastpm_store_get_lagrangian_position(FastPMStore * p, ptrdiff_t index, double pos[3])
{
    pos[0] = p->q[index][0];
    pos[1] = p->q[index][1];
    pos[2] = p->q[index][2];
}

double fastpm_store_get_mass(FastPMStore * p, ptrdiff_t index)
{
    /* total mass is the sum of the base and the extra */
    if(p->mass) {
        return p->meta.M0 + p->mass[index];
    } else {
        return p->meta.M0;
    }
}

static ptrdiff_t
_alignsize(ptrdiff_t size)
{
    /* align sizes */
    return ((size + 1024) / 1024) * 1024;
}

void
fastpm_store_init_details(FastPMStore * p,
                  const char * name,
                  size_t np_upper,
                  FastPMColumnTags attributes,
                  enum FastPMMemoryLocation loc,
                  const char * file,
                  const int line)
{
    p->mem = _libfastpm_get_gmem();

    /* only set the name if name is not NULL; this is to allow setting the name before calling init.
     * */
    if(name) {
        strcpy(p->name, name);
    }

    p->attributes = attributes;

    p->np = 0;
    p->np_upper = np_upper;

    /* clear the column pointers. */
    memset(p->columns, 0, sizeof(p->columns));
    memset(p->_column_info, 0, sizeof(p->_column_info));

    /* clear the meta */
    memset(&p->meta, 0, sizeof(p->meta));
    int it;
    p->_base = NULL;

    /* first loop initialize the _column_info struct for corresponding items in the pointer union;
     * second loop initialize the pointers */
    #define DEFINE_COLUMN(column, attr_, dtype_, nmemb_) \
        { \
            int ci = FASTPM_STORE_COLUMN_INDEX(column); \
            if(attr_ != (1 << ci)) { fastpm_raise(-1, "attr and column are out of order for %s\n", # column ); } \
            strcpy(p->_column_info[ci].dtype, dtype_);  \
            strcpy(p->_column_info[ci].name, # column); \
            p->_column_info[ci].elsize = sizeof(p->column[0]);  \
            p->_column_info[ci].nmemb = nmemb_;  \
            p->_column_info[ci].membsize = sizeof(p->column[0]) / nmemb_;  \
            p->_column_info[ci].attribute = attr_;  \
            p->_column_info[ci].pack = pack_any;  \
            p->_column_info[ci].unpack = unpack_any;  \
            p->_column_info[ci].from_double = NULL;  \
            p->_column_info[ci].to_double = NULL;  \
        }

    #define COLUMN_INFO(column) (p->_column_info[FASTPM_STORE_COLUMN_INDEX(column)])

    DEFINE_COLUMN(x, COLUMN_POS, "f8", 3);
    DEFINE_COLUMN(q, COLUMN_Q, "f4", 3);
    DEFINE_COLUMN(v, COLUMN_VEL, "f4", 3);
    DEFINE_COLUMN(acc, COLUMN_ACC, "f4", 3);
    DEFINE_COLUMN(dx1, COLUMN_DX1, "f4", 3);
    DEFINE_COLUMN(dx2, COLUMN_DX2, "f4", 3);
    DEFINE_COLUMN(dv1, COLUMN_DV1, "f4", 3);
    DEFINE_COLUMN(aemit, COLUMN_AEMIT, "f4", 1);
    DEFINE_COLUMN(rho, COLUMN_DENSITY, "f4", 1);
    DEFINE_COLUMN(potential, COLUMN_POTENTIAL, "f4", 1);
    DEFINE_COLUMN(tidal, COLUMN_TIDAL, "f4", 6);
    DEFINE_COLUMN(id, COLUMN_ID, "i8", 1);
    DEFINE_COLUMN(pgdc, COLUMN_PGDC, "f4", 3);
    DEFINE_COLUMN(mask, COLUMN_MASK, "i1", 1);
    DEFINE_COLUMN(minid, COLUMN_MINID, "i8", 1);
    DEFINE_COLUMN(task, COLUMN_TASK, "i4", 1);
    DEFINE_COLUMN(length, COLUMN_LENGTH, "i4", 1);
    DEFINE_COLUMN(rdisp, COLUMN_RDISP, "f4", 6);
    DEFINE_COLUMN(vdisp, COLUMN_VDISP, "f4", 6);
    DEFINE_COLUMN(rvdisp, COLUMN_RVDISP, "f4", 9);
    DEFINE_COLUMN(mass, COLUMN_MASS, "f4", 1);

    COLUMN_INFO(x).to_double = to_double_f8;
    COLUMN_INFO(v).to_double = to_double_f4;
    COLUMN_INFO(rho).to_double = to_double_f4;
    COLUMN_INFO(dx1).to_double = to_double_f4;
    COLUMN_INFO(dx2).to_double = to_double_f4;
    COLUMN_INFO(dv1).to_double = to_double_f4;
    COLUMN_INFO(acc).to_double = to_double_f4;
    COLUMN_INFO(mass).to_double = to_double_f4;

    COLUMN_INFO(rho).from_double = from_double_f4;
    COLUMN_INFO(acc).from_double = from_double_f4;
    COLUMN_INFO(pgdc).from_double = from_double_f4;
    COLUMN_INFO(dx1).from_double = from_double_f4;
    COLUMN_INFO(dx2).from_double = from_double_f4;
    COLUMN_INFO(dv1).from_double = from_double_f4;
    COLUMN_INFO(potential).from_double = from_double_f4;
    COLUMN_INFO(tidal).from_double = from_double_f4;

    ptrdiff_t size = 0;
    ptrdiff_t offset = 0;
    for(it = 0; it < 2; it ++) {
        int ci;
        for(ci = 0; ci < 32; ci ++) {
            if(it == 0) {
                size += ((attributes & p->_column_info[ci].attribute) != 0) * (_alignsize(p->_column_info[ci].elsize * np_upper)); \
            } else { \
                if(attributes & p->_column_info[ci].attribute) { \
                    p->columns[ci] = (void*) (((char*) p->_base) + offset); \
                    offset += _alignsize(p->_column_info[ci].elsize * np_upper); \
                } else { \
                    p->columns[ci] = NULL; \
                } \
            }

        }
        if(it == 0) {
            p->_base = fastpm_memory_alloc_details(p->mem, "FastPMStore", size, loc, file, line);
            /* zero out all memory */
            memset(p->_base, 0, size);
        }
    };
}

void
fastpm_store_set_name(FastPMStore * p, const char * name)
{
    strcpy(p->name, name);
}

size_t
fastpm_store_init_evenly_details(FastPMStore * p,
    const char * name, size_t np_total, FastPMColumnTags attributes, double alloc_factor, MPI_Comm comm,
    const char * file,
    const int line)
{
    /* allocate for np_total cross all */
    /* name means name of species*/
    int NTask;
    MPI_Comm_size(comm, &NTask);

    size_t np_upper = (size_t)(1.0 * np_total / NTask * alloc_factor);

    MPI_Bcast(&np_upper, 1, MPI_LONG, 0, comm);
    fastpm_store_init_details(p, name, np_upper, attributes, FASTPM_MEMORY_HEAP, file, line);
    return 0;
}

size_t
fastpm_store_get_np_total(FastPMStore * p, MPI_Comm comm)
{
    long long np = p->np;
    MPI_Allreduce(MPI_IN_PLACE, &np, 1, MPI_LONG_LONG, MPI_SUM, comm);
    return np;
}

size_t
fastpm_store_get_mask_sum(FastPMStore * p, MPI_Comm comm)
{
    long long np = 0;
    ptrdiff_t i;
    for(i = 0; i < p->np; i ++) {
        np += p->mask[i] != 0;
    }
    MPI_Allreduce(MPI_IN_PLACE, &np, 1, MPI_LONG_LONG, MPI_SUM, comm);
    return np;
}

void
fastpm_packing_plan_init(FastPMPackingPlan * plan, FastPMStore * p, FastPMColumnTags attributes)
{
    int ci;
    int i = 0;
    plan->elsize = 0;
    plan->attributes = attributes;
    for (ci = 0; ci < 32; ci ++) {
        if (!(p->_column_info[ci].attribute & attributes)) continue;
        plan->_ci[i] = ci;
        plan->_offsets[ci] = plan->elsize;
        ptrdiff_t elsize = p->_column_info[ci].elsize;
        plan->elsize += elsize;
        plan->_column_info[ci] = p->_column_info[ci];
        i++;
    }
    /* Pad the elsize to 8 bytes. This ensures anything goes to the MPI wire is 8 byte aligned,
     * and minimizes the chances of hitting an implementation bug. */
    if (plan->elsize % 8 != 0) {
        plan->elsize += (8 - plan->elsize % 8);
    }
    plan->Ncolumns = i;
}

void
fastpm_packing_plan_pack(FastPMPackingPlan * plan,
            FastPMStore * p, ptrdiff_t i, void * packed)
{
    int t;
    memset(packed, 0, plan->elsize);
    for (t = 0; t < plan->Ncolumns; t ++) {
        int ci = plan->_ci[t];
        ptrdiff_t offset = plan->_offsets[ci];
        plan->_column_info[ci].pack(p, i, ci,
            ((char*) packed) + offset);
    }
}

void
fastpm_packing_plan_unpack(FastPMPackingPlan * plan,
            FastPMStore * p, ptrdiff_t i, void * packed)
{
    int t;
    for (t = 0; t < plan->Ncolumns; t ++) {
        int ci = plan->_ci[t];
        fastpm_packing_plan_unpack_ci(plan, ci, p, i, packed);
    }
}

/* unpack a single column from the offset in packed data. */
void
fastpm_packing_plan_unpack_ci(FastPMPackingPlan * plan, int ci,
            FastPMStore * p, ptrdiff_t i, void * packed)
{
    ptrdiff_t offset = plan->_offsets[ci];
    plan->_column_info[ci].unpack(p, i, ci,
        ((char*) packed) + offset);
}

int
fastpm_store_find_column_id(FastPMStore * p, FastPMColumnTags attribute)
{
    int ci;
    for (ci = 0; ci < 32; ci ++) {
        if (p->_column_info[ci].attribute == attribute) {
            return ci;
        }
    }
    fastpm_raise(-1, "Unknown column %x", attribute);
    return -1;
}

void
fastpm_store_destroy(FastPMStore * p)
{
    fastpm_memory_free(p->mem, p->_base);
}

static void permute(void * data, int np, size_t elsize, int * ind) {
    void * tmp = malloc(elsize * np);
    if(!tmp) {
        fastpm_raise(-1, "No memory for permuting\n");
    }
    int i;
    for(i = 0; i < np; i ++) {
        memcpy(((char*) tmp) + i * elsize, ((char*) data) + ind[i] * elsize, elsize);
    }
    memcpy(data, tmp, np * elsize);
    free(tmp);
}

void fastpm_store_permute(FastPMStore * p, int * ind)
{
    int c;
    for(c = 0; c < 32; c ++) {
        if(!p->columns[c]) continue;
        permute(p->columns[c], p->np, p->_column_info[c].elsize, ind);
    }
}

static FastPMStore *  _fastpm_store_sort_store;
static int (*_fastpm_store_sort_cmp_func)(const int i1, const int i2, FastPMStore * p);

int _sort_by_id_cmpfunc(const void * p1, const void * p2)
{
    const int * i1 = (const int*) p1;
    const int * i2 = (const int*) p2;

    return _fastpm_store_sort_cmp_func(*i1, *i2, _fastpm_store_sort_store);
}

int
FastPMLocalSortByID(const int i1,
                    const int i2,
                    FastPMStore * p)
{
    int v1 = (p->id[i1] < p->id[i2]);
    int v2 = (p->id[i1] > p->id[i2]);

    return v2 - v1;
}

/* sort a store locally with in the MPI rank.
 *
 * cmp_func is called with three arguments.
 * */
void
fastpm_store_sort(FastPMStore * p,
        int (*cmp_func)(const int i1, const int i2, FastPMStore * p))
{
    int * arg = fastpm_memory_alloc(p->mem, "Temp", sizeof(int) * p->np, FASTPM_MEMORY_HEAP);
    int i;
    for(i = 0; i < p->np; i ++) {
        arg[i] = i;
    }
    /* FIXME: copy some version of qsort_r */
    _fastpm_store_sort_store = p;
    _fastpm_store_sort_cmp_func = cmp_func;
    qsort(arg, p->np, sizeof(arg[0]), _sort_by_id_cmpfunc);
    fastpm_store_permute(p, arg);
    fastpm_memory_free(p->mem, arg);
}

void
fastpm_store_wrap(FastPMStore * p, double BoxSize[3])
{
    int i;
    int d;
    for(i = 0; i < p->np; i ++) {
        for(d = 0; d < 3; d ++) {
            double n = abs(p->x[i][d] / BoxSize[d]);

            double x1 = remainder(p->x[i][d], BoxSize[d]);

            while(x1 < 0) x1 += BoxSize[d];
            while(x1 > BoxSize[d]) x1 -= BoxSize[d];
            p->x[i][d] = x1;

            if(n > 10000) {
                double q[3];
                if(fastpm_store_has_q(p)) {
                    fastpm_store_get_q_from_id(p, p->id[i], q);
                }
                fastpm_raise(-1, "Particle at %g %g %g (q = %g %g %g) is too far from the bounds. Wrapping failed.\n",
                        p->x[i][0],
                        p->x[i][1],
                        p->x[i][2],
                        q[0], q[1], q[2]
                );
            }
        }
    }
}

int
FastPMTargetPM (FastPMStore * p, ptrdiff_t i, PM * pm)
{
    double pos[3];
    fastpm_store_get_position(p, i, pos);
    return pm_pos_to_rank(pm, pos);
}

int
fastpm_store_decompose(FastPMStore * p,
    fastpm_store_target_func target_func,
    void * data, MPI_Comm comm)
{
    if(fastpm_store_get_np_total(p, comm) == 0) return 0 ;

    VALGRIND_CHECK_MEM_IS_DEFINED(p->x, sizeof(p->x[0]) * p->np);

    FastPMPackingPlan plan[1];

    fastpm_packing_plan_init(plan, p, p->attributes);

    size_t elsize = plan->elsize;

    size_t Nsend_limit = 1000 * 1024 * 1024 / elsize;     // this large number effectively prevents throttling

    int NTask, ThisTask;

    MPI_Comm_rank(comm, &ThisTask);
    MPI_Comm_size(comm, &NTask);

    if(p->np > p->np_upper) {
        fastpm_raise(-1, "Particle buffer overrun detected np = %td > np_upper %td.\n", p->np, p->np_upper);
    }

    /* do a bincount; offset by -1 because -1 is for self */
    int * count = calloc(NTask + 1, sizeof(int));
    int * offsets = calloc(NTask + 1, sizeof(int));
    int * sendcount = count + 1;
    int * recvcount = malloc(sizeof(int) * (NTask));
    int * recvoffset = malloc(sizeof(int) * (NTask));
    int * sendoffset = malloc(sizeof(int) * (NTask));

    int incomplete = 1;
    int iter = 0;
    /* terminate based on incomplete for throttling*/
    while(1) {
        incomplete = 0;
        int * target = fastpm_memory_alloc(p->mem, "Target", sizeof(int) * p->np, FASTPM_MEMORY_HEAP);

        size_t Nsend_all = 0;
        ptrdiff_t i;
        for(i = 0; i < p->np; i ++) {
            target[i] = target_func(p, i, data);
            if(ThisTask == target[i]) {
                target[i] = -1;
            } else {
                Nsend_all++;
                /* Throttling: never send more than this many particles. */
                if(Nsend_all >= Nsend_limit) {
                    incomplete = 1;
                    target[i] = -1;
                } else {
                }
            }
        }

        for(i = 0; i < NTask + 1; i ++) {
            count[i] = 0;
            offsets[i] = 0;
        }
        for(i = 0; i < NTask; i ++) {
            sendoffset[i] = 0;
            recvoffset[i] = 0;
        }

        for(i = 0; i < p->np; i ++) {
            sendcount[target[i]] ++;
        }
        cumsum(offsets, count, NTask + 1);

        int * arg = fastpm_memory_alloc(p->mem, "PermArg", sizeof(int) * p->np, FASTPM_MEMORY_HEAP);
        for(i = 0; i < p->np; i ++) {
            int offset = offsets[target[i] + 1] ++;
            arg[offset] = i;
        }

        fastpm_store_permute(p, arg);

        VALGRIND_CHECK_MEM_IS_DEFINED(p->x, sizeof(p->x[0]) * p->np);

        fastpm_memory_free(p->mem, arg);
        fastpm_memory_free(p->mem, target);

        MPI_Alltoall(sendcount, 1, MPI_INT,
                     recvcount, 1, MPI_INT,
                     comm);

        size_t Nsend = cumsum(sendoffset, sendcount, NTask);
        size_t Nrecv = cumsum(recvoffset, recvcount, NTask);

        volatile size_t neededsize = p->np + Nrecv - Nsend;

        if(neededsize > p->np_upper) {
            fastpm_ilog(INFO, "Need %td particles on rank %d; %td allocated\n", neededsize, ThisTask, p->np_upper);
        }

        if(MPIU_Any(comm, neededsize > p->np_upper)) {
            goto fail_oom;
        }

        p->np -= Nsend;

        void * send_buffer = fastpm_memory_alloc(p->mem, "SendBuf", elsize * Nsend, FASTPM_MEMORY_HEAP);
        void * recv_buffer = fastpm_memory_alloc(p->mem, "RecvBuf", elsize * Nrecv, FASTPM_MEMORY_HEAP);

        {
            double nmin, nmax, nmean, nstd;

            MPIU_stats(comm, elsize * Nsend, "<>-s", &nmin, &nmax, &nmean, &nstd);
            fastpm_info("Send buffer size : min=%g max=%g mean=%g, std=%g bytes", nmin, nmax, nmean, nstd);
            MPIU_stats(comm, elsize * Nrecv, "<>-s", &nmin, &nmax, &nmean, &nstd);
            fastpm_info("Recv buffer size : min=%g max=%g mean=%g, std=%g bytes", nmin, nmax, nmean, nstd);
        }

        int ipar = 0;
        int j;
        for(i = 0; i < NTask; i ++) {
            for(j = sendoffset[i]; j < sendoffset[i] + sendcount[i]; j ++, ipar++) {
                fastpm_packing_plan_pack(plan, p, j + p->np, (char*) send_buffer + ipar * elsize);
            }
        }

        size_t Nsendsum;
        size_t Nsendallsum;
        MPI_Allreduce(&Nsend, &Nsendsum, 1, MPI_LONG, MPI_SUM, comm);
        MPI_Allreduce(&Nsend_all, &Nsendallsum, 1, MPI_LONG, MPI_SUM, comm);
        fastpm_info("Decomposition iter %d,  exchange of %td particles; need %td", iter, Nsendsum, Nsendallsum);

        MPI_Datatype PTYPE;
        MPI_Type_contiguous(elsize, MPI_BYTE, &PTYPE);
        MPI_Type_commit(&PTYPE);

        MPI_Alltoallv_sparse(
                send_buffer, sendcount, sendoffset, PTYPE,
                recv_buffer, recvcount, recvoffset, PTYPE,
                comm);

        MPI_Type_free(&PTYPE);

        ipar = 0;
        for(i = 0; i < NTask; i ++) {
            for(j = recvoffset[i]; j < recvoffset[i] + recvcount[i]; j++, ipar ++) {
                fastpm_packing_plan_unpack(plan, p, j + p->np, (char*) recv_buffer + ipar * elsize);
            }
        }

        fastpm_memory_free(p->mem, recv_buffer);
        fastpm_memory_free(p->mem, send_buffer);

        p->np += Nrecv;
        iter++;
        if(MPIU_Any(comm, incomplete)) continue;
        else break;
    }
    free(recvcount);
    free(recvoffset);
    free(sendoffset);
    free(offsets);
    free(count);

    return 0;

    fail_oom:
        free(recvcount);
        free(recvoffset);
        free(sendoffset);
        free(offsets);
        free(count);
    return -1;
}

int fastpm_store_has_q(FastPMStore * p)
{
    return p->meta._q_size != 0;
}

void
fastpm_store_get_q_from_id(FastPMStore * p, uint64_t id, double q[3])
{
    ptrdiff_t pabs[3];

    int d;
    id = id % p->meta._q_size;
    for(d = 0; d < 3; d++) {
        pabs[d] = id / p->meta._q_strides[d];
        id -= pabs[d] * p->meta._q_strides[d];
    }

    for(d = 0; d < 3; d++) {
        q[d] = pabs[d] * p->meta._q_scale[d];
        q[d] += p->meta._q_shift[d];
    }
}//if all 0 drop

void
fastpm_store_get_iq_from_id(FastPMStore * p, uint64_t id, ptrdiff_t pabs[3])
{
    /* integer version of the initial position q
    i.e. the lattice coordinates of the cells. */
    int d;
    for(d = 0; d < 3; d++) {
        pabs[d] = id / p->meta._q_strides[d];
        id -= pabs[d] * p->meta._q_strides[d];
    }
}

void
fastpm_store_fill(FastPMStore * p, PM * pm, double * shift, ptrdiff_t * Nc)
{
    /* fill p with a uniform grid, respecting domain given by pm. use a subsample ratio.
     * (every subsample grid points) */
    if(Nc == NULL) {
        Nc = pm_nmesh(pm);
    }
    int d;
    p->np = 1;
    for(d = 0; d < 3; d++) {
        int start = pm->IRegion.start[d] * Nc[d] / pm->Nmesh[d];
        int end = (pm->IRegion.start[d] + pm->IRegion.size[d]) * Nc[d] / pm->Nmesh[d];
        p->np *= end - start;
    }
    if(p->np > p->np_upper) {
        fastpm_raise(-1, "Need %td particles; %td allocated\n", p->np, p->np_upper);
    }
    ptrdiff_t ptr = 0;

    for(d = 0; d < 3; d ++) {
        if(shift)
            p->meta._q_shift[d] = shift[d];
        else
            p->meta._q_shift[d] = 0;

        p->meta._q_scale[d] = pm->BoxSize[d] / Nc[d];
    }

    p->meta._q_size = Nc[0] * Nc[1] * Nc[2];
    p->meta._q_strides[0] = Nc[1] * Nc[2];
    p->meta._q_strides[1] = Nc[2];
    p->meta._q_strides[2] = 1;

    PMXIter iter;
    for(pm_xiter_init(pm, &iter);
       !pm_xiter_stop(&iter);
        pm_xiter_next(&iter)){
        ptrdiff_t pabs_start[3];
        ptrdiff_t pabs_end[3];
        ptrdiff_t ii, jj, kk;

        for(d = 0; d < 3; d ++) {
            pabs_start[d] = iter.iabs[d] * Nc[d] / pm->Nmesh[d];
            pabs_end[d] = (iter.iabs[d] + 1) * Nc[d] / pm->Nmesh[d];
        }
        for(ii = pabs_start[0]; ii < pabs_end[0]; ii ++)
        for(jj = pabs_start[1]; jj < pabs_end[1]; jj ++)
        for(kk = pabs_start[2]; kk < pabs_end[2]; kk ++) {
            ptrdiff_t pabs[3] = {ii, jj, kk};

            uint64_t id = pabs[2] * p->meta._q_strides[2] +
                          pabs[1] * p->meta._q_strides[1] +
                          pabs[0] * p->meta._q_strides[0] ;

            if(p->id) p->id[ptr] = id;
            if(p->mask) p->mask[ptr] = 0;

            fastpm_store_get_q_from_id(p, id, &p->x[ptr][0]);

            if(p->q) {
                /* set q if it is allocated. */
                for(d = 0; d < 3; d ++) {
                    p->q[ptr][d] = p->x[ptr][d];
                }
            }
            ptr ++;
        }
    }
    if(ptr != p->np) {
        fastpm_raise(-1, "This is an internal error, particle number mismatched with grid. %td != %td, allocsize=%td, shape=(%td %td %td)\n",
            ptr, p->np, pm->allocsize,
            pm->IRegion.size[0],
            pm->IRegion.size[1],
            pm->IRegion.size[2]
            );
    }
    p->meta.a_x = p->meta.a_v = 0.;
}

void
fastpm_store_summary(FastPMStore * p,
        FastPMColumnTags attribute,
        MPI_Comm comm,
        const char * fmt,
        ...)
{
    va_list va;
    va_start(va, fmt);

    int ci = fastpm_store_find_column_id(p, attribute);
    size_t nmemb = p->_column_info[ci].nmemb;

    double rmin[nmemb], rmax[nmemb], rsum1[nmemb], rsum2[nmemb];

    if (NULL == p->_column_info[ci].to_double) {
        fastpm_raise(-1, "Column %s didnot set to_double virtual function\n",
            p->_column_info[ci].name);
    }
    int d;
    for(d = 0; d < nmemb; d ++) {
        rmin[d] = 1e20;
        rmax[d] = -1e20;
        rsum1[d] = 0;
        rsum2[d] = 0;
    }
    ptrdiff_t i;

#pragma omp parallel
    {
        double tmin[nmemb], tmax[nmemb], tsum1[nmemb], tsum2[nmemb];
        int d;
        for(d = 0; d < nmemb; d ++) {
            tmin[d] = 1e20;
            tmax[d] = -1e20;
            tsum1[d] = 0;
            tsum2[d] = 0;
        }
        #pragma omp for
        for(i = 0; i < p->np; i ++) {
            int d;
            for(d =0; d < nmemb; d++) {
                double value = p->_column_info[ci].to_double(p, i, ci, d);
                tsum1[d] += value;
                tsum2[d] += value * value;
                tmin[d] = fmin(tmin[d], value);
                tmax[d] = fmax(tmax[d], value);
            }
        }
        #pragma omp critical
        {
            int d;
            for(d =0; d < nmemb; d++) {
                rsum1[d] += tsum1[d];
                rsum2[d] += tsum2[d];
                rmin[d] = fmin(rmin[d], tmin[d]);
                rmax[d] = fmax(rmax[d], tmax[d]);
            }
        }
    }
    uint64_t Ntot = p->np;

    MPI_Allreduce(MPI_IN_PLACE, rsum1, nmemb, MPI_DOUBLE, MPI_SUM, comm);
    MPI_Allreduce(MPI_IN_PLACE, rsum2, nmemb, MPI_DOUBLE, MPI_SUM, comm);
    MPI_Allreduce(MPI_IN_PLACE, rmin, nmemb, MPI_DOUBLE, MPI_MIN, comm);
    MPI_Allreduce(MPI_IN_PLACE, rmax, nmemb, MPI_DOUBLE, MPI_MAX, comm);
    MPI_Allreduce(MPI_IN_PLACE, &Ntot,   1, MPI_LONG,  MPI_SUM, comm);

    for(i = 0; i < strlen(fmt); i ++) {
        void * r = va_arg(va, void *);
        double * dr = (double * ) r;

        for(d = 0; d < 3; d++) {
            switch(fmt[i]) {
                case '-':
                    dr[d] = rsum1[d] / Ntot;
                    break;
                case '<':
                    dr[d] = rmin[d];
                    break;
                case '>':
                    dr[d] = rmax[d];
                    break;
                case 's':
                    dr[d] = sqrt(rsum2[d] / Ntot - pow(rsum1[d] / Ntot, 2));
                    break;
                case 'S':
                    dr[d] = sqrt(1.0 * Ntot / (Ntot - 1.)) * sqrt(rsum2[d] / Ntot - pow(rsum1[d] / Ntot, 2));
                    break;
                case 'v':
                    dr[d] = (rsum2[d] / Ntot - pow(rsum1[d] / Ntot, 2));
                    break;
                case 'V':
                    dr[d] = (1.0 * Ntot / (Ntot - 1.)) * (rsum2[d] / Ntot - pow(rsum1[d] / Ntot, 2));
                    break;
                default:
                    fastpm_raise(-1, "Unknown format str. Use '<->sSvV'\n");
            }
        }
    }
    va_end(va);
}

void
fastpm_store_steal(FastPMStore * p, FastPMStore * po, FastPMColumnTags attributes)
{
    int c;
    for(c = 0; c < 32; c ++) {
        if (!(p->_column_info[c].attribute & attributes)) continue;
        po->columns[c] = p->columns[c];
    }

    po->np = p->np;
    po->meta = p->meta;
}

static void
_fastpm_store_copy(FastPMStore * p, ptrdiff_t start, FastPMStore * po, ptrdiff_t offset, size_t ncopy)
{
    if(ncopy + start > p->np) {
        fastpm_raise(-1, "Copy out of bounds from source FastPMStore: asking for %td but has %td\n", ncopy + start, p->np);
    }
    if(ncopy + offset > po->np_upper) {
        fastpm_raise(-1, "Not enough storage in target FastPMStore: asking for %td but has %td\n", ncopy + offset, po->np_upper);
    }

    int c;
    for(c = 0; c < 32; c ++) {
        if(!po->columns[c]) continue;
        size_t elsize = po->_column_info[c].elsize;

        memcpy(po->columns[c] + offset * elsize,
                p->columns[c] + start * elsize, elsize * ncopy);

    }
    po->np = offset + ncopy;
    po->meta = p->meta;
}

void
fastpm_store_copy(FastPMStore * p, FastPMStore * po)
{
    _fastpm_store_copy(p, 0, po, 0, p->np);
}

void
fastpm_store_take(FastPMStore * p, ptrdiff_t i, FastPMStore * po, ptrdiff_t j)
{
    _fastpm_store_copy(p, i, po, j, 1);
}

/* extends p by extra. */
void
fastpm_store_extend(FastPMStore * p, FastPMStore * extra)
{
    _fastpm_store_copy(extra, 0, p, p->np, extra->np);
}

void
fastpm_store_fill_subsample_mask(FastPMStore * p,
        double fraction,
        FastPMParticleMaskType * mask,
        MPI_Comm comm)
{
    gsl_rng * random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
    int ThisTask;
    MPI_Comm_rank(comm, &ThisTask);

    /* set uncorrelated seeds */
    double seed=1231584; //FIXME: set this properly.
    gsl_rng_set(random_generator, seed);
    int d;

    for(d = 0; d < ThisTask * 8; d++) {
        seed = 0x7fffffff * gsl_rng_uniform(random_generator);
    }

    gsl_rng_set(random_generator, seed);

    memset(mask, 0, p->np * sizeof(mask[0]));

    ptrdiff_t i;
    for(i=0; i < p->np; i++) {
        double rand_i = gsl_rng_uniform(random_generator);
        int flag = fraction > 1 || rand_i <= fraction;
        mask[i] = flag;
    }

    gsl_rng_free(random_generator);
}

void
fastpm_store_fill_subsample_mask_every_dim(FastPMStore * p,
                                              int every, /* take 1 every 'every' per dimension */
                                              FastPMParticleMaskType * mask)
{
    if(!fastpm_store_has_q(p)) {
        /* This can be relaxed by using a less strict subsample algorithm, e.g. subsample by a hash of ID. */
        fastpm_raise(-1, "Subsample is not supported if the store does not have meta.q.");
    }
    /* UNUSED */
    memset(mask, 0, p->np * sizeof(mask[0]));

    ptrdiff_t i, d;
    for(i = 0; i < p->np; i++) {
        ptrdiff_t pabs[3];   //pabs is the iq index. move out of loop?
        uint64_t id = p->id[i];
        fastpm_store_get_iq_from_id(p, id, pabs);

        int flag = 1;
        for(d = 0; d < 3; d++) {
            flag *= !(pabs[d] % every);
        }
        mask[i] = flag;
    }
}

/*
 * Create a subsample, keeping only those with mask == True;
 *
 * if po is NULL, only return number of items.
 * */
size_t
fastpm_store_subsample(FastPMStore * p, FastPMParticleMaskType * mask, FastPMStore * po)
{
    ptrdiff_t i;
    ptrdiff_t j;

    j = 0;
    for(i = 0; i < p->np; i ++) {
        if(!mask[i]) continue;
        /* just counting */
        if(po == NULL) { j++; continue; }
        /* avoid memcpy of same address if we are doing subsample inplace */
        if(p == po && j == i) {j ++; continue; }

        int c;
        for(c = 0; c < 32; c ++) {
            if (!po->columns[c]) continue;

            size_t elsize = po->_column_info[c].elsize;
            memcpy(po->columns[c] + j * elsize, p->columns[c] + i * elsize, elsize);
        }

        j ++;
    }

    if(po) {
        po->np = j;
        po->meta = p->meta;   ///????
    }
    return j;
}

#if 0
static ptrdiff_t
binary_search(double foo, double a[], size_t n) {
    ptrdiff_t left = 0, right = n;
    ptrdiff_t mid;
    /* find the right side, because hist would count the number < edge, not <= edge */
    if(a[left] > foo) {
        return 0;
    }
    if(a[right - 1] <= foo) {
        return n;
    }
    while(right - left > 1) {
        mid = left + ((right - left - 1) >> 1);
        double pivot = a[mid];
        if(pivot > foo) {
            right = mid + 1;
            /* a[right - 1] > foo; */
        } else {
            left = mid + 1;
            /* a[left] <= foo; */
        }
    }
    return left;
}
#endif
/* this is cumulative */
void
fastpm_store_histogram_aemit_sorted(FastPMStore * store,
        int64_t * hist,
        double * aedges,
        size_t nedges,
        MPI_Comm comm)
{
    ptrdiff_t i;

    int64_t * hist1 = malloc(sizeof(hist1[0]) * (nedges + 1));

    memset(hist1, 0, sizeof(hist1[0]) * (nedges + 1));

#pragma omp parallel
    {
        /* FIXME: use standard reduction with OpenMP 4.7 */

        int64_t * hist2 = malloc(sizeof(hist2[0]) * (nedges + 1));
        memset(hist2, 0, sizeof(hist2[0]) * (nedges + 1));

        int iedge = 0;

        /* this works because openmp will send each thread at most 1
         * chunk with the static scheduling; thus a thread never sees
         * out of order aemit */
        #pragma omp for schedule(static)
        for(i = 0; i < store->np; i ++) {
            while(iedge < nedges && store->aemit[i] >= aedges[iedge]) iedge ++;

//           int ibin = binary_search(store->aemit[i], aedges, nedges);
//           if(ibin != iedge) abort();

            hist2[iedge] ++;
        }

        #pragma omp critical
        {
            for(i = 0; i < nedges + 1; i ++) {
                hist1[i] += hist2[i];
            }
        }

        free(hist2);
    }

    MPI_Allreduce(MPI_IN_PLACE, hist1, nedges + 1, MPI_LONG, MPI_SUM, comm);

    for(i = 0; i < nedges + 1; i ++) {
        hist[i] += hist1[i];
    }

    free(hist1);
}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner