text:
/*

   Copyright 2016  Amazon.com, Inc. or its affiliates. All Rights Reserved.

   Licensed under the Apache License, Version 2.0 (the "License"). You may not use this file except in compliance with the License. A copy of the License is located at

   http://aws.amazon.com/apache2.0/

   or in the "license" file accompanying this file. This file 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 NNTYPES_H
#include <vector>
#include <set>

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <algorithm>
#include <netcdf>
#ifndef __NVCC__
#include <tuple>
#include <json/json.h>
#endif
#include <sys/time.h>
#include <cmath>

class NNDataSetBase;
class NNLayer;
class NNNetwork;
class NNWeight;

// Activates step by step CPU validation
#define VALIDATION
#ifdef VALIDATION
extern "C"
{
    #include <cblas.h>
}
#endif

static const float NN_VERSION       = 0.81f;
static const float MIN_ERROR        = 1.0e-12f;
static const float MIN_ACTIVATION   = 0.000001f;
static const float MAX_ACTIVATION   = 0.999999f;
static const float MAX_VALUE        = 999999999999999.0f;

template <typename T> struct GpuBuffer;

enum
{
    DefaultBatch    = 512
};

enum TrainingMode
{
    SGD = 0,
    Momentum = 1,
    AdaGrad = 2,
    Nesterov = 3,
    RMSProp = 4,
    AdaDelta = 5,
};

ostream& operator<< (ostream& out, const TrainingMode& e);

enum ErrorFunction
{
    L1,
    L2,
    CrossEntropy,
    ScaledMarginalCrossEntropy,
};

ostream& operator<< (ostream& out, const ErrorFunction& e);

enum Activation {
    Sigmoid,
    Tanh,
    RectifiedLinear,
    Linear,
    ParametricRectifiedLinear,
    SoftPlus,
    SoftSign,
    SoftMax,
    ReluMax,
    LinearMax,
};

ostream& operator<< (ostream& out, const Activation& a);

enum WeightInitialization
{
    Xavier,
    CaffeXavier,
    Gaussian,
    Uniform,
    UnitBall,
    Constant
};

ostream& operator<< (ostream& out, const WeightInitialization& w);

enum PoolingFunction {
    None,
    Max,
    Average,
    Stochastic,
    LocalContrastNormalization,
    LocalResponseNormalization,
    GlobalTemporal,
};

ostream& operator<< (ostream& out, const PoolingFunction& p);

#include "kernels.h"
#include "GpuSort.h"
#include "NNEnum.h"
#include "NNWeight.h"
#include "NNLayer.h"
#include "NNNetwork.h"

int MPI_Bcast_string(string& s);

struct NNDataSetDimensions
{
    uint32_t _dimensions;
    uint32_t _width;
    uint32_t _height;
    uint32_t _length;
};

struct NNDataSetBase {

    string                   _name;                          // Dataset name
    NNDataSetEnums::DataType _dataType;                      // Dataset type (see above enum)
    uint32_t                 _attributes;                    // Dataset characteristics (see NNDataSetEnum::Attributes in NNEnum.h)
    uint32_t                 _examples;                      // Number of examples
    uint32_t                 _dimensions;                    // Dimensionality of data set
    uint32_t                 _width;                         // Dataset x dimension
    uint32_t                 _height;                        // Dataset y dimension
    uint32_t                 _length;                        // Dataset z dimension
    uint32_t                 _stride;                        // Stride between examples
    NNDataSetEnums::Sharding _sharding;                      // Sharding of dataset for parallel execution
    uint32_t                 _minX;                          // Beginning of local X sharding for model parallel execution
    uint32_t                 _maxX;                          // End of local X sharding for model parallel execution
    uint64_t                 _sparseDataSize;                // Total sparse datapoints
    uint32_t                 _maxSparseDatapoints;           // Maximum observed sparse datapoints per example
    NNFloat                  _sparseDensity;                 // Overall sparse density (0.0 - 1.0)
    vector<uint64_t>         _vSparseStart;                  // Vector of sparse datapoint starts per example
    GpuBuffer<uint64_t>*     _pbSparseStart;                 // GPU copy of _vSparseStart
    vector<uint64_t>         _vSparseEnd;                    // Vector of sparse datapoint ends per example
    GpuBuffer<uint64_t>*     _pbSparseEnd;                   // GPU copy of _vSparseEnd
    vector<uint32_t>         _vSparseIndex;                  // Vector of sparse indices
    GpuBuffer<uint32_t>*     _pbSparseIndex;                 // GPU copy of _vSparseIndex
    GpuBuffer<NNFloat>*      _pbDenoisingRandom;             // Denoising randoms

    // Transposed sparse lookup for sparse backpropagation
    vector<uint64_t>         _vSparseDatapointCount;
    vector<uint32_t>         _vSparseTransposedStart;
    uint32_t                 _sparseTransposedIndices;
    GpuBuffer<uint32_t>*     _pbSparseTransposedStart;
    GpuBuffer<uint32_t>*     _pbSparseTransposedEnd;
    GpuBuffer<uint32_t>*     _pbSparseTransposedIndex;

    // States
    bool                     _bDenoising;
    bool                     _bDirty;
    uint32_t                 _batch;

    NNDataSetBase();
    NNDataSetDimensions GetDimensions();
    uint32_t GetExamples() { return _examples; };

    virtual bool SaveNetCDF(const string& fname) = 0;
    virtual bool WriteNetCDF(netCDF::NcFile& nfc, const string& fname, const uint32_t n) = 0;
    virtual ~NNDataSetBase() = 0;
    virtual void RefreshState(uint32_t batch) = 0;
    virtual bool Shard(NNDataSetEnums::Sharding sharding) = 0;
    virtual bool UnShard() = 0;
    virtual vector<tuple<uint64_t, uint64_t> > getMemoryUsage() = 0;
    virtual bool CalculateSparseDatapointCounts() = 0;
    virtual bool GenerateSparseTransposedMatrix(uint32_t batch, NNLayer* pLayer) = 0;
    virtual bool CalculateSparseTransposedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer) = 0;
    virtual bool CalculateSparseTransposedDenoisedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer) = 0;
    virtual bool CalculateSparseTransposedWeightGradient(NNFloat alpha, NNFloat beta, uint32_t m, uint32_t n, NNFloat* pDelta, NNFloat* pWeightGradient) = 0;
    virtual bool SetDenoising(bool flag) = 0;
    virtual bool GenerateDenoisingData() = 0;
    virtual bool LoadInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual bool LoadSparseInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual bool LoadSparseDenoisedInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual bool CalculateSparseZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta = (NNFloat)0.0) = 0;
    virtual bool CalculateSparseDenoisedZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta = (NNFloat)0.0) = 0;
    virtual float CalculateL1Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual float CalculateL2Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual float CalculateCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual float CalculateScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual float CalculateMultinomialCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual float CalculateMultinomialScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit) = 0;
    virtual bool CalculateL1OutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta) = 0;
    virtual bool CalculateCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta) = 0;
    virtual bool CalculateScaledMarginalCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta) = 0;
    virtual bool CalculateOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta) = 0;
};

ostream& operator<< (ostream& out, NNDataSetEnums::Attributes& a);
ostream& operator<< (ostream& out, NNDataSetEnums::Kind& k);
ostream& operator<< (ostream& out, NNDataSetEnums::DataType& t);
ostream& operator<< (ostream& out, NNDataSetEnums::Sharding& s);

template<typename T> class NNDataSet : public NNDataSetBase {
public:
    friend class NNetwork;
    friend class NNLayer;
    friend vector<NNDataSetBase*> LoadNetCDF(const string& fname);
    friend bool SaveNetCDF(const string& fname, vector<NNDataSetBase*> vDataSet);

private:

    vector<T>               _vData;
    GpuBuffer<T>*           _pbData;
    vector<T>               _vSparseData;
    GpuBuffer<T>*           _pbSparseData;
    GpuBuffer<T>*           _pbSparseTransposedData;

    // Force constructor private
    NNDataSet(const string& fname, uint32_t n);
    bool Rename(const string& name);
    bool SaveNetCDF(const string& fname);
    bool WriteNetCDF(netCDF::NcFile& nfc, const string& fname, const uint32_t n);
    void RefreshState(uint32_t batch) {}
    bool Shard(NNDataSetEnums::Sharding sharding);
    bool UnShard();
    vector<tuple<uint64_t, uint64_t> > getMemoryUsage();
    bool CalculateSparseDatapointCounts();
    bool GenerateSparseTransposedMatrix(uint32_t batch, NNLayer* pLayer);
    bool CalculateSparseTransposedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer);
    bool CalculateSparseTransposedDenoisedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer);
    bool CalculateSparseTransposedWeightGradient(NNFloat alpha, NNFloat beta, uint32_t m, uint32_t n, NNFloat* pDelta, NNFloat* pWeightGradient);
    bool SetDenoising(bool flag);
    bool GenerateDenoisingData();
    bool LoadInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    bool LoadSparseInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    bool LoadSparseDenoisedInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    bool CalculateSparseZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta);
    bool CalculateSparseDenoisedZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta);
    float CalculateL1Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    float CalculateL2Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    float CalculateCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    float CalculateScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    float CalculateMultinomialCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    float CalculateMultinomialScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit);
    bool CalculateL1OutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta);
    bool CalculateCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta);
    bool CalculateScaledMarginalCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta);
    bool CalculateOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta);

public:

    ~NNDataSet();
    void Shuffle();
    T GetDataPoint(uint32_t n, uint32_t x, uint32_t y = 0, uint32_t z = 0);
    bool SetDataPoint(T v, uint32_t n, uint32_t x, uint32_t y = 0, uint32_t z = 0);
    uint32_t GetSparseDataPoints(uint32_t n);
    uint32_t GetSparseIndex(uint32_t n, uint32_t i);
    bool SetSparseIndex(uint32_t n, uint32_t i, uint32_t v);
    T GetSparseDataPoint(uint32_t n, uint32_t i);
    bool SetSparseDataPoint(uint32_t n, uint32_t i, T v);
};

template<typename T> bool NNDataSet<T>::LoadInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    kLoadInputUnit(position, batch, stride, pUnit, _pbData->_pDevData);
    return true;
}

template<typename T> bool NNDataSet<T>::LoadSparseInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Boolean)
        kLoadSparseInputUnit(position, batch, stride, pUnit, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData);
    else
        kLoadSparseAnalogInputUnit(position, batch, stride, pUnit, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData);
    return true;
}

template<typename T> bool NNDataSet<T>::LoadSparseDenoisedInputUnit(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Boolean)
        kLoadSparseDenoisedInputUnit(position, batch, stride, pUnit, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbDenoisingRandom->_pDevData);
    else
        kLoadSparseAnalogDenoisedInputUnit(position, batch, stride, pUnit, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, _pbDenoisingRandom->_pDevData);
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateSparseZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta)
{
    if (_attributes & NNDataSetEnums::Boolean)
        kCalculateSparseZ(position, batch, stride, pWeight, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, pUnit, beta);
    else
        kCalculateSparseAnalogZ(position, batch, stride, pWeight, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, pUnit, beta);
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateSparseDenoisedZ(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pWeight, NNFloat* pUnit, NNFloat beta)
{
    if (_attributes & NNDataSetEnums::Boolean)
        kCalculateSparseDenoisedZ(position, batch, stride, pWeight, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbDenoisingRandom->_pDevData, pUnit, beta);
    else
        kCalculateSparseAnalogDenoisedZ(position, batch, stride, pWeight, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, _pbDenoisingRandom->_pDevData, pUnit, beta);
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateSparseTransposedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer)
{
    // Rebuild sparse data table if dataset changed
    if (_bDirty || (batch != _batch))
    {
        GenerateSparseTransposedMatrix(batch, pLayer);
    }

    // Initialize transposed sparse offsets
    _pbSparseTransposedEnd->Copy(_pbSparseTransposedStart->_pDevData);

    // Call appropriate matrix generation kernel
    if (_attributes & NNDataSetEnums::Boolean)
        kCalculateSparseTransposedMatrix(position, batch, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData);
    else
        kCalculateSparseTransposedAnalogMatrix(position, batch, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData, _pbSparseTransposedData->_pDevData);

    return true;
}

template<typename T> bool NNDataSet<T>::CalculateSparseTransposedDenoisedMatrix(uint32_t position, uint32_t batch, NNLayer* pLayer)
{

    // Rebuild sparse data table if dataset changed
    if (_bDirty || (batch != _batch))
    {
        GenerateSparseTransposedMatrix(batch, pLayer);
    }

    // Initialize transposed sparse offsets
    _pbSparseTransposedEnd->Copy(_pbSparseTransposedStart->_pDevData);

    // Call appropriate matrix generation kernel
    if (_attributes & NNDataSetEnums::Boolean)
        kCalculateSparseTransposedDenoisedMatrix(position, batch, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbDenoisingRandom->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData);
    else
        kCalculateSparseTransposedAnalogDenoisedMatrix(position, batch, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, _pbDenoisingRandom->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData, _pbSparseTransposedData->_pDevData);

#if 0
    vector<uint32_t> vSparseTransposedStart(53120);
    vector<uint32_t> vSparseTransposedEnd(53120);
    _pbSparseTransposedStart->Download(&vSparseTransposedStart[0]);
    _pbSparseTransposedEnd->Download(&vSparseTransposedEnd[0]);
    for (uint32_t i = 0; i < 53120; i++)
    printf("%6u %9u %9u %9u %9u\n", i, vSparseTransposedStart[i], vSparseTransposedEnd[i], vSparseTransposedEnd[i] - vSparseTransposedStart[i], (uint32_t)_vSparseDatapointCount[i]);
    exit(-1);
#endif
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateSparseTransposedWeightGradient(NNFloat alpha, NNFloat beta, uint32_t m, uint32_t n, NNFloat* pDelta, NNFloat* pWeightGradient)
{
    if (_attributes & NNDataSetEnums::Boolean)
        kCalculateSparseTransposedWeightGradient(alpha, beta, m, n, _pbSparseTransposedStart->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData, pDelta, pWeightGradient);
    else
        kCalculateSparseTransposedAnalogWeightGradient(alpha, beta, m, n, _pbSparseTransposedStart->_pDevData, _pbSparseTransposedEnd->_pDevData, _pbSparseTransposedIndex->_pDevData, _pbSparseTransposedData->_pDevData, pDelta, pWeightGradient);
    return true;
}

template<typename T> float NNDataSet<T>::CalculateL1Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        if (_attributes & NNDataSetEnums::Boolean)
           return kCalculateSparseL1Error(position, batch, stride, pUnit,
                  _pbSparseStart->_pDevData,
                  _pbSparseEnd->_pDevData,
                  _pbSparseIndex->_pDevData,
                  bSparseIgnoreZero);
        else
           return kCalculateSparseAnalogL1Error(position, batch, stride, pUnit,
                  _pbSparseStart->_pDevData,
                  _pbSparseEnd->_pDevData,
                  _pbSparseIndex->_pDevData,
                  _pbSparseData->_pDevData,
                  bSparseIgnoreZero);
    }
    else
        return kCalculateL1Error(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> float NNDataSet<T>::CalculateL2Error(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        if (_attributes & NNDataSetEnums::Boolean)
            return kCalculateSparseL2Error(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData,
                   bSparseIgnoreZero);
        else
            return kCalculateSparseAnalogL2Error(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData,
                   _pbSparseData->_pDevData,
                   bSparseIgnoreZero);
    }
    else
        return kCalculateL2Error(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> float NNDataSet<T>::CalculateCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        return kCalculateSparseCrossEntropyError(position, batch, stride, pUnit,
               _pbSparseStart->_pDevData,
               _pbSparseEnd->_pDevData,
               _pbSparseIndex->_pDevData,
               bSparseIgnoreZero);
    }
    else
        return kCalculateCrossEntropyError(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> float NNDataSet<T>::CalculateScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        return kCalculateSparseScaledMarginalCrossEntropyError(position, batch, stride, pUnit,
               _pbSparseStart->_pDevData,
               _pbSparseEnd->_pDevData,
               _pbSparseIndex->_pDevData,
               bSparseIgnoreZero);
    }
    else
        return kCalculateScaledMarginalCrossEntropyError(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> float NNDataSet<T>::CalculateMultinomialCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        if (_attributes & NNDataSetEnums::Boolean)
        {
            return kCalculateSparseMultinomialCrossEntropyError(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData);
        }
        else
            return kCalculateSparseAnalogMultinomialCrossEntropyError(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData,
                   _pbSparseData->_pDevData);
    }
    else
        return kCalculateMultinomialCrossEntropyError(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> float NNDataSet<T>::CalculateMultinomialScaledMarginalCrossEntropyError(uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        if (_attributes & NNDataSetEnums::Boolean)
            return kCalculateSparseMultinomialScaledMarginalCrossEntropyError(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData);
        else
            return kCalculateSparseAnalogMultinomialScaledMarginalCrossEntropyError(position, batch, stride, pUnit,
                   _pbSparseStart->_pDevData,
                   _pbSparseEnd->_pDevData,
                   _pbSparseIndex->_pDevData,
                   _pbSparseData->_pDevData);
    }
    else
        return kCalculateMultinomialScaledMarginalCrossEntropyError(position, batch, stride, pUnit, _pbData->_pDevData);
}

template<typename T> bool NNDataSet<T>::CalculateL1OutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        kCalculateSparseL1OutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, bSparseIgnoreZero);
    }
    else
        kCalculateL1OutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbData->_pDevData);
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        kCalculateSparseCrossEntropyOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, bSparseIgnoreZero);
    }
    else
        kCalculateCrossEntropyOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbData->_pDevData);
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateScaledMarginalCrossEntropyOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta)
{
    if (_attributes & NNDataSetEnums::Sparse)
    {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        kCalculateSparseScaledMarginalCrossEntropyOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, bSparseIgnoreZero);
    }
    else
    {
        kCalculateScaledMarginalCrossEntropyOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbData->_pDevData);
    }
    return true;
}

template<typename T> bool NNDataSet<T>::CalculateOutputDelta(Activation activation, uint32_t position, uint32_t batch, uint32_t stride, NNFloat* pUnit, NNFloat* pDelta)
{
    if (_attributes & NNDataSetEnums::Sparse) {
        bool bSparseIgnoreZero = _attributes & NNDataSetEnums::SparseIgnoreZero;
        if (_attributes & NNDataSetEnums::Boolean)
        {
            kCalculateSparseOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, bSparseIgnoreZero);
        }
        else
        {
            kCalculateSparseAnalogOutputDelta(activation, position, batch, stride, pUnit,  pDelta, _pbSparseStart->_pDevData, _pbSparseEnd->_pDevData, _pbSparseIndex->_pDevData, _pbSparseData->_pDevData, bSparseIgnoreZero);
        }
    }
    else
    {
        kCalculateOutputDelta(activation, position, batch, stride, pUnit, pDelta, _pbData->_pDevData);
    }
    return true;
}

vector<NNDataSetBase*> LoadNetCDF(const string& fname);
bool SaveNetCDF(const string& fname, vector<NNDataSetBase*> vDataset);
vector<NNDataSetBase*> LoadImageData(const string& fname);
vector<NNDataSetBase*> LoadCSVData(const string& fname);
vector<NNDataSetBase*> LoadJSONData(const string& fname);
vector<NNDataSetBase*> LoadAudioData(const string& name);

#define NNTYPES_H
#endif


text:
/* block/gsl_block_char.h
 *
 * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_BLOCK_CHAR_H__
#define __GSL_BLOCK_CHAR_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_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

struct gsl_block_char_struct
{
  size_t size;
  char *data;
};

typedef struct gsl_block_char_struct gsl_block_char;

GSL_EXPORT gsl_block_char *gsl_block_char_alloc (const size_t n);
GSL_EXPORT gsl_block_char *gsl_block_char_calloc (const size_t n);
GSL_EXPORT void gsl_block_char_free (gsl_block_char * b);

GSL_EXPORT int gsl_block_char_fread (FILE * stream, gsl_block_char * b);
GSL_EXPORT int gsl_block_char_fwrite (FILE * stream, const gsl_block_char * b);
GSL_EXPORT int gsl_block_char_fscanf (FILE * stream, gsl_block_char * b);
GSL_EXPORT int gsl_block_char_fprintf (FILE * stream, const gsl_block_char * b, const char *format);

GSL_EXPORT int gsl_block_char_raw_fread (FILE * stream, char * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_char_raw_fwrite (FILE * stream, const char * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_char_raw_fscanf (FILE * stream, char * b, const size_t n, const size_t stride);
GSL_EXPORT int gsl_block_char_raw_fprintf (FILE * stream, const char * b, const size_t n, const size_t stride, const char *format);

GSL_EXPORT size_t gsl_block_char_size (const gsl_block_char * b);
GSL_EXPORT char * gsl_block_char_data (const gsl_block_char * b);

__END_DECLS

#endif /* __GSL_BLOCK_CHAR_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


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

#include "optimization.h"

/*****************************************************************************/
/* vector routines */

double vec_dot_prod(const double *v1, const double *v2, const int N)
{
    int i;
    double nn = 0;
    for (i=0;i<N;i++)
    {
        nn += v1[i] * v2[i];
    }
    return nn;
}

double vec_norm(const double *v1, const int N)
{
    return sqrt(vec_dot_prod(v1, v1, N));
}

/*****************************************************************************/
/* Wolfe conditions */

int check_armijo(double new_misfit, double old_misfit,
                 double dir_deriv, double step_len)
{
    /* returns 1 if Armijo condition satisfied, 0 otherwise */
    if (new_misfit <= old_misfit + WOLFE_C1 * step_len * dir_deriv)
    {
        /* step length okay */
        return 1;
    }
    return 0;
}

int check_curvature(double dir_deriv_old, double dir_deriv_new)
{
    /* returns 1 if curvature satisfied, 0 otherwise */
    if (dir_deriv_new >= WOLFE_C2 * dir_deriv_old)
    {
        return 1;
    }
    return 0;
}

int check_strong_curvature(double dir_deriv_old, double dir_deriv_new)
{
    /* returns 1 if strong curvature satisfied, 0 otherwise */
    if (dir_deriv_new >=  WOLFE_C2 * dir_deriv_old &&
        dir_deriv_new <= -WOLFE_C2 * dir_deriv_old)
    {
        return 1;
    }
    return 0;
}

/*****************************************************************************/

double initial_step_length(search_dir sd, double alpha_old, double dir_deriv_old, double dir_deriv_new)
{
    /* calculate initial step length for a line search step         */
    /* inputs:                                                      */
    /*  sd: search direction, CG, STEEPEST, or LBFGS                */
    /*  alpha_old: previous (accepted) step length                  */
    /*  dir_deriv_old: directional derivative at previous step      */
    /*      (dot product of previous gradient and search direction) */
    /*  dir_deriv_new: directional derivative at current step       */
    /*      (dot product of current gradient and search direction)  */

    if (sd == CG || sd == STEEPEST)
    {
        /* see Nocedal  & Wright 2006, p. 59 */
        return alpha_old * (dir_deriv_old / dir_deriv_new);
    }

    /* quasi-Newton methods are well scaled use steplength of 1.0 for L-BFGS */
    return 1.0;

}

/*****************************************************************************/
/* steepest descent direction */

void steepest_descent_direction(const double *grad, double *p, const int N)
{
    int i;
    for (i=0;i<N;i++)
    {
        p[i] = -grad[i];
    }
}

/*****************************************************************************/
/* conjugate gradient direction */

void conjugate_grad_direction(double *grad_old, double *grad_new, double *p_old, double *p_new, int N)
{
    /*
    Determine a search direction by the conjugate gradient method

    Uses PR-FR conjugate gradient method to calculate a conjugate direection,
    see Nocedal & Wright 2006, p. 123.

    All vectors should be double precision of length N.
    grad_old: previous gradient
    grad_new: current gradient
    p_old: previous search direction
    p_new: new search direction (output, must be allocated already)
    */

    /* calculate beta values for both PR and FR methods */
    double beta_lower = 0;
    double beta_upper_FR = 0;
    double beta_upper_PR = 0;

    int i;
    for (i=0;i<N;i++)
    {
        beta_lower += grad_old[i] * grad_old[i];
        beta_upper_FR += grad_new[i] * grad_new[i];
        beta_upper_PR += grad_new[i] * (grad_new[i] - grad_old[i]);
    }

    double beta_FR = beta_upper_FR / beta_lower;
    double beta_PR = beta_upper_PR / beta_lower;

    /* select beta based on criteria given in Nocedal and Wright 2006 */
    double beta = 0;

    if (beta_PR < -beta_FR)
    {
        beta = -beta_FR;
    }
    else if (fabs(beta_PR) <= beta_FR)
    {
        beta = beta_PR;
    }
    else if (beta_PR > beta_FR)
    {
        beta = beta_FR;
    }

    /* calculate new search direction */
    for (i=0;i<N;i++)
    {
        p_new[i] = -grad_new[i] + beta * p_old[i];
    }
}

/*****************************************************************************/
/* L-BFGS direction */

void LBFGS_direction(double *grad_new, double *grad_old, double *model_new, double *model_old, double *p_new, int N, int m, int k)
{
    /* call with the last parameter negative to clean up local variables */

    int ii, jj;

    static double *s = NULL;
    static double *y = NULL;

    static double *q = NULL;
    static double *aa = NULL;
    static double *rho = NULL;
    static int *ind = NULL;

    double *r = p_new;

    /* allocate memory */
    /* note: opt_allocate is the same as malloc, but checks to make sure    */
    /* memory allocate was successful                                       */
    if (!s)
    {
        s = opt_allocate(sizeof(*s) * N * m);
    }

    if (!y)
    {
        y = opt_allocate(sizeof(*y) * N * m);
    }

    if (!q)
    {
        q = opt_allocate(sizeof(*q) * N);
    }

    if (!aa)
    {
        aa = opt_allocate(sizeof(*aa) * m);
    }

    if (!rho)
    {
        rho = opt_allocate(sizeof(*rho) * m);
    }

    if (!ind)
    {
        ind = opt_allocate(sizeof(*ind) * m);
        for (ii=0;ii<m;ii++)
        {
            ind[ii] = ii * N;
        }
    }

    /* j is the min of k and m */
    int j = (k < m) ? k : m;

    /* position to write new y and s vectors */
    int pos = ((k-1) % m) * N;

    /* write new y and s vectors */
    for (ii=0;ii<N;ii++)
    {
        y[ii+pos] = grad_new[ii] - grad_old[ii];
        s[ii+pos] = model_new[ii] - model_old[ii];
    }

    /* calculate rho */
    rho[pos / N] = 1.0 / vec_dot_prod(y+pos, s+pos, N);

    /* calculate gamma */
    double gamma = vec_dot_prod(s+pos, y+pos, N) / vec_dot_prod(y+pos, y+pos, N);

    if (k<m)
    {
        /* shift ind vector so ind[0] is still the index
        of the oldest element */
        int first_ele = ind[0];
        for (ii=0;ii<m-1;ii++)
        {
            ind[ii] = ind[ii+1];
        }
        ind[m-1] = first_ele;
    }

    /* assign initial q value */
    for (jj=0;jj<N;jj++)
    {
        q[jj] = grad_new[jj];
    }

    /*begin backwards loop */
    for (ii=j-1;ii>=0;ii--)
    {
        aa[ii] = rho[ind[ii] / N] * vec_dot_prod(s+ind[ii], q, N);

        for (jj=0;jj<N;jj++)
        {
            q[jj] = q[jj] - aa[ii] * (y+ind[ii])[jj];
        }
    }

    /* assign r */
    for (jj=0;jj<N;jj++)
    {
        r[jj] = gamma * q[jj];
    }

    /* begin forwards loop */
    double beta;
    for (ii=0;ii<j;ii++)
    {
        beta = rho[ind[ii] / N] * vec_dot_prod(y+ind[ii], r, N);

        for (jj=0;jj<N;jj++)
        {
            r[jj] = r[jj] + (s+ind[ii])[jj] * (aa[ii] - beta);
        }
    }

    /* negate r */
    for (jj=0;jj<N;jj++)
    {
        r[jj] = -r[jj];
    }

}

/*****************************************************************************/
#include <lapacke.h>
void newton_direction(const double *gradient, const double *hessian, const int N, double *p_new)
{
    /********************************************************/
    /* not implemented yet - use steepest descent for now   */
    //fprintf(stderr,"Newton direction not implemented yet\n");
    //fprintf(stderr,"Reverting to steepest descent\n");
    //steepest_descent_direction(gradient, p_new, N);
    /********************************************************/

    /* make a temp copy of the gradient*/
    double *temp_g = opt_allocate(sizeof(*temp_g) * N);
    double *temp_h = opt_allocate(sizeof(*temp_h) * N * N);

    memcpy(temp_g, gradient, sizeof(*temp_h) * N );

    /* need to check if Hessian is PD and modify if not                 */
    /* Note that this repeated performs a Cholesky factorization on     */
    /* the Hessian and will be expensive if N is large                  */
    /* In those cases, a modified Cholesky factorization (see Nocedal   */
    /* and Wright 2006, pp. 52-54) may be a better approach             */
    double beta = sqrt(DBL_EPSILON);
    int pd = hessian_modification(hessian, temp_h, N, beta);

    /* if we found a pd Hessian approximation, solve system for p       */
    /* otherwise, return steepest descent direction                     */
    if (pd != 1)
    {
        fprintf(stderr,"Could not find a positive definite "
                       "Hessian approximation\n");
        fprintf(stderr,"Returning steepest descent direction\n");
    }
    else
    {
        /* solve system for p: Hp = -g */
        int nrhs = 1;
        int lda = N;
        int ldb = N;
        LAPACKE_dpotrs(LAPACK_COL_MAJOR, 'L', N, nrhs, temp_h, lda, temp_g, ldb);
    }

    /* negate gradient or solution to system */
    int i;
    for (i = 0;i<N;i++)
        p_new[i]  = -temp_g[i];

    free(temp_h);
    free(temp_g);
}

#define HESSMOD_MAX_IT 100
int hessian_modification(const double *hessian, double *upper_factor, int N, double beta)
{
    /* beta may need to be adjusted based on the scaling of the problem */

    /* find minimum diagonal element of the hessian */
    int i;
    double min_diag = hessian[0];
    for (i=0;i<N;i++)
    {
        double diag = hessian[i*(N+1)];
        min_diag = (diag < min_diag) ? diag : min_diag;
    }

    /* find initial tau */
    double tau;
    if (min_diag > 0)
    {
        /* if it is already pd, min_diag will be positive */
        tau = 0.0;
    }
    else
    {
        tau = - min_diag + beta;
    }

    /* add multiples of the identity until hessian is positive definite */
    int k;
    int pd;
    int lda = N;
    for (k=0;k<HESSMOD_MAX_IT;k++)
    {
        /* copy hessian and add identity multiple */
        /* upper_factor will change every loop because of the dpotrf call */
        memcpy(upper_factor, hessian, sizeof(*hessian) * N *N);
        for (i=0;i<N;i++)
        {
            upper_factor[i * (N+1)] = hessian[i * (N+1)] + tau;
        }

        /* try cholesky factorization */
        pd = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', N, upper_factor, lda);
        if (pd == 0)
        {
            /* positive definite */
            return 1;
        }
        /* increment tau */
        tau = 2.0 * tau > beta ? 2.0 * tau : beta;

    }
    /* return 0 if not positive definite */
    return 0;

}

/*****************************************************************************/
/* function for performing line search */

void line_search(double *x_in,
                       double *x_out,
                       const int N,
                       optimization_parameters *opt_par,
                       void *aux
                       )
{

    double (*misfit)(const double *, const int, void *) = opt_par->objective;
    void (*grad)(const double *, double *, int, void *) = opt_par->grad;
    void (*hess)(const double *, double *, int, void *) = opt_par->hess;

    int success = 0;
    int i,j;

    /* note: opt_allocate is the same as malloc, but checks to make sure    */
    /* memory allocate was successful                                       */

    double *gradient = opt_allocate(sizeof(*gradient) * N);
    double *grad_old = opt_allocate(sizeof(*grad_old) * N);

    double *p = opt_allocate(sizeof(*p) * N);
    double *p_old = opt_allocate(sizeof(*p_old) * N);

    double *x_old = opt_allocate(sizeof(*x_old) * N);
    double *x_trial = opt_allocate(sizeof(*x_trial) * N);

    double *hessian = NULL;
    if (opt_par->sd == NEWTON)
        hessian = opt_allocate(sizeof(*hessian) * N * N);

    /* initial alpha_old and dir_deriv_old to keep the compiler from complaining,   */
    /* but these initialization values will never be used                           */
    double alpha_old = 1;
    double alpha_new;

    double dir_deriv_old = 1;
    double dir_deriv_new;

    /* copy x_in to x_old */
    for (i=0;i<N;i++)
    {
        x_out[i] = x_in[i];
    }

    /* evaluate function and gradient at initial point */
    double f_old = misfit(x_out, N, aux);
    grad(x_out, gradient, N, aux);

    /* check grad norm for exit condition */
    if (vec_norm(gradient, N) < opt_par->stopping_tolerance)
    {
        success = 1;
        free(p);
        free(gradient);
        return;
    }

    int m = opt_par->LBFGS_mem;

    /* start main loop */
    search_dir sd;
    int max_it = opt_par->max_iterations;
    for (j=0;j<max_it;j++)
    {
        printf("Starting line search, iteration %d\n",j);
        if (j == 0 && opt_par->sd != NEWTON)
        {
            /* first iteration - use steepest descent */
            /* get steepest descent direction */
            steepest_descent_direction(gradient, p, N);
            sd = STEEPEST;
            dir_deriv_new = vec_dot_prod(gradient, p, N);
            alpha_new = opt_par->initial_step_length;
        }
        else
        {
            if (opt_par->sd == LBFGS)
            {
                LBFGS_direction(gradient, grad_old, x_out, x_old, p, N, m, j);

                sd = LBFGS;
                dir_deriv_new = vec_dot_prod(gradient, p, N);
            }
            else if (opt_par->sd == CG)
            {
                conjugate_grad_direction(grad_old, gradient, p_old, p, N);
                sd = CG;
                dir_deriv_new = vec_dot_prod(gradient, p, N);
            }
            else if (opt_par->sd == NEWTON)
            {
                /* compute hessian */
                hess(x_out, hessian, N, aux);
                newton_direction(gradient, hessian, N, p);
                sd = NEWTON;
                dir_deriv_new = vec_dot_prod(gradient, p, N);
            }

            if (opt_par->sd == STEEPEST || (opt_par->angle_restart_allowed && check_angle_restart(gradient, p, N)))
            {
                steepest_descent_direction(gradient, p, N);
                sd = STEEPEST;
                dir_deriv_new = vec_dot_prod(gradient, p, N);
            }

            alpha_new = initial_step_length(sd, alpha_old, dir_deriv_old, dir_deriv_new);
        }

        if (opt_par->bracket_only || j==0)
        {
            alpha_new = ls_bracket(x_out, p, N, opt_par, alpha_new, aux);
        }
        else
        {
            alpha_new = back_track(x_out, p, N, opt_par, alpha_new, aux);
        }

        for (i=0;i<N;i++)
        {
            p_old[i] = p[i];
            grad_old[i] = gradient[i];
            x_old[i] = x_out[i];
        }

        if (!opt_par->backtrack_only && !opt_par->bracket_only && j!=0)
        {
            /* update model */
            for (i=0;i<N;i++)
            {
                x_trial[i] = x_out[i] + alpha_new * p[i];
            }

            /* evaluate gradient at new point */
            grad(x_trial, gradient, N, aux);

            /* check (strong) curvature */
            double ppi = vec_dot_prod(gradient, p, N);
            if (!check_strong_curvature(dir_deriv_new, ppi))
            {
                //alpha_new = initial_step_length(sd, alpha_old, dir_deriv_old, dir_deriv_new);
                alpha_new = ls_bracket(x_out, p, N, opt_par, alpha_new, aux);
            }

        }

        alpha_old = alpha_new;
        dir_deriv_old = dir_deriv_new;

        /* calculate new x value */
        for (i=0; i<N; i++)
        {
            x_out[i] = x_out[i] + alpha_new * p[i];
        }

        /* evaluate gradient at new point */
        grad(x_out, gradient, N, aux);

        /* check grad norm for exit condition */
        if (vec_norm(gradient, N) < opt_par->stopping_tolerance)
        {
            success = 1;
            break;
        }

    }

    printf("No. iterations: %d\n",j);

    free(p);
    free(p_old);
    free(gradient);
    free(grad_old);
    free(x_trial);
    free(x_old);

    if (hessian)
        free(hessian);

    if (!success)
    {
        fprintf(stderr,"Error: line search failed to find a solution\n");
        fprintf(stderr,"Point returned is not stationary\n");
    }

}

/*****************************************************************************/

double ls_zoom(const double *x,
               const double *p,
               const int N,
               optimization_parameters *opt_par,
               double a_lo,
               double a_hi,
               void *aux
               )
{
    double (*misfit)(const double *, const int, void *) = opt_par->objective;
    void (*grad)(const double *, double *, int, void *) = opt_par->grad;

    int success = 0;

    double bb_star = 0.0;

    double mf;
    double *gradient = opt_allocate(sizeof(*gradient) * N);
    double *x_update = opt_allocate(sizeof(*x_update) * N);

    /* calculate initial misfit and gradient */
    mf = misfit(x, N, aux);
    grad(x, gradient, N, aux);

    double phi0 = mf;
    double phi_prime0 = vec_dot_prod(gradient, p, N);

    double phi_aj;
    double phi_prime_aj;

    double phi_a_lo;

    double mean_step;

    double tol = 0;

    /* start iterations */
    double a_j;
    int i, j;
    for (i=0;i<opt_par->bracket_zoom_iterations;i++)
    {
        printf("Bracket successful, starting zoom, iteration %d\n",i);
        printf("a_lo, a_hi: %f, %f\n", a_lo, a_hi);
        fflush(stdout);

        /* check if a_lo and a_hi are too close together */
        mean_step = 0.5 * (a_lo + a_hi);
        if (fabs(a_lo - a_hi) / mean_step < 0.01)
        {
            bb_star = mean_step;
            success = 1;
            break;
        }
        if (mean_step < 1e-6) /* <-- this should really depend on the scaling of the problem */
        {
            success = 0;
            break;
        }

        a_j = ls_cubic_interpolation(x, p, N, opt_par, a_lo, a_hi, aux);
        //a_j = (a_lo + a_hi) / 2.0;

        tol = mean_step * 0.01; /* one percent of the mean step */
        //fprintf(stderr, "Trial step: %f\n",a_j);
        if (fabs(a_j - a_lo) < tol || fabs(a_j - a_hi) < tol)
        {
            /* too close to boundary of interval */
            //fprintf(stderr,"ls_zoom: too close to boundary\n");
            //success = 0;
            //break;
            a_j = mean_step;
        }

        if (isnan(a_j) || isinf(a_j))
        {
            fprintf(stderr,"nan/inf value\n");
            success = 0;
            break;
        }

        /* evaluate function at a_lo (don't need dir. derivative) */
        for (j=0;j<N;j++)
        {
            x_update[j] = x[j] + a_lo * p[j];
        }
        phi_a_lo = misfit(x_update, N, aux);

        /* evaluate function at a_j. directional derivative will be         */
        /* calculated later if needed                                       */
        for (j=0;j<N;j++)
        {
            x_update[j] = x[j] + a_j * p[j];
        }
        phi_aj = misfit(x_update, N, aux);

        printf("Trial step length %E misfit: %E\n", a_j, phi_aj);

        if (phi_aj > phi0 + WOLFE_C1 * a_j * phi_prime0
            || phi_aj >= phi_a_lo)
        {
            /* trial step length violates Armijo or                         */
            /* trial step evaluate to a function value greater than a_lo    */
            /* need to reduce step length                                   */
            a_hi = a_j;
        }
        else
        {
            /* trial step length satisfies Armijo and function value is     */
            /* less than a_lo                                               */
            /* check curvature                                              */
            grad(x_update, gradient, N, aux);
            phi_prime_aj = vec_dot_prod(gradient, p, N);
            printf("Trial step length %E dir deriv: %E\n", a_j, phi_prime_aj);
            if (fabs(phi_prime_aj) <= -WOLFE_C2 * phi_prime0)
            {
                /* trial step length satisfies strong curvature condition   */
                /* accept trial step length                                 */
                bb_star = a_j;
                success = 1;
                break;
            }

            if (phi_prime_aj * (a_hi - a_lo) >= 0)
            {
                a_hi = a_lo;
            }

            /* satisfies Armijo, but not curvature                          */
            /* increase step length                                         */
            a_lo = a_j;

        }
    }

    free(x_update);
    free(gradient);

    if (!success)
    {
        fprintf(stderr,"Line search zoom failed to return an "
        "appropriate step length\n");
        //exit(-90);
        return -1;
    }

    return bb_star;
}

/*****************************************************************************/

double ls_cubic_interpolation(const double *x,
               const double *p,
               const int N,
               optimization_parameters *opt_par,
               const double bound1,
               const double bound2,
               void *aux
               )
{
    /* calculates a new trial step length from cubic interpolation  */
    /* interpolation is based on two function values and the        */
    /* corresponding directional derivatives                        */
    /*                                                              */
    /* This interpolation requires two function evaluations and     */
    /* two gradient evaluations. If the function or gradient is     */
    /* expensive to evaluate, there are probably better ways to     */
    /* calculate an appropriate step length                         */
    /*                                                              */
    /* see Nocedal and Wright 2006, equation 3.59                   */
    /*                                                              */
    /* inputs:                                                      */
    /*      x: current point in parameter space                     */
    /*      p: search direction                                     */
    /*      N: number of dimensions                                 */
    /*      misfit: function that takes point and returns objective */
    /*              function value. x and N as defined above        */
    /*      grad: function that takes point and returns gradient    */
    /*            value (returned as parameter)                     */
    /*      bound1, bound2: bounds of step length search            */
    /*                                                              */
    /*  new trial step length is returned                           */

    double (*misfit)(const double *, const int, void *) = opt_par->objective;
    void (*grad)(const double *, double *, int, void *) = opt_par->grad;

    double *gradient = opt_allocate(sizeof(*gradient) * N);
    double *x_update = opt_allocate(sizeof(*x_update) * N);

    double f1, f2, g1, g2;

    int i;

    /* evaluate objective function and directional derivative at bound1 */
    for (i=0;i<N;i++)
    {
        x_update[i] = x[i] + bound1 * p[i];
    }
    f1 = misfit(x_update, N, aux);
    grad(x_update, gradient, N, aux);
    g1 = vec_dot_prod(gradient, p, N);

    printf("misfit at    %E: %E\n", bound1, f1);
    printf("dir deriv at %E: %E\n", bound1, g1);

    /* evaluate objective function and directional derivative at bound2 */
    for (i=0;i<N;i++)
    {
        x_update[i] = x[i] + bound2 * p[i];
    }
    f2 = misfit(x_update, N, aux);
    grad(x_update, gradient, N, aux);
    g2 = vec_dot_prod(gradient, p, N);

    printf("misfit at    %E: %E\n", bound2, f2);
    printf("dir deriv at %E: %E\n", bound2, g2);

    /*calculate cubic approximateion */
    double d1 = g1 + g2 - 3 * ((f1 - f2)/(bound1 - bound2));

    /* note that both bound1 and bound2 should always be positive   */
    /* they should never be equal                                   */
    /* d2_sign is sign(bound2 - bound1)                             */
    double d2_sign = bound2 > bound1 ? 1.0 : -1.0;
    if ((bound2 - bound1) > 0)
        d2_sign = 1.0;
    else if ((bound2 - bound1) < 0)
        d2_sign = -1.0;
    else
        d2_sign = 0.0;

    double d2 = d2_sign * sqrt(d1 * d1 - g1 * g2);

    double a_new = bound2 - (bound2 - bound1) * ((g2 + d2 - d1)/(g2 - g1 + 2 * d2));

    /* the minimum is always either at a_new as calculated above orj    */
    /* at one of the end points. Check those here.                      */
    //a_new = (a_new < bound1) ? a_new : bound1;
    //a_new = (a_new < bound2) ? a_new : bound2;

    free(x_update);
    free(gradient);
    return a_new;
}

/*****************************************************************************/

double ls_bracket(const double *x,
               const double *p,
               const int N,
               optimization_parameters *opt_par,
               double alpha0,
               void *aux
               )
{

    double (*misfit)(const double *, const int, void *) = opt_par->objective;
    void (*grad)(const double *, double *, int, void *) = opt_par->grad;

    int success = 0;
    double b_star;

    int steplength_max = 50;
    double beta_max = 5e30;
    double beta0 = 0;
    double beta1 = alpha0; /* initial step length */

    double *gradient = opt_allocate(sizeof(*gradient) * N);
    double *x_update = opt_allocate(sizeof(*x_update) * N);

    double p0;      /* phi(0)  */
    double pp0;     /* phi'(0) */
    /* evaluate objective and dir. derivative at initial point */
    p0 = misfit(x, N, aux);
    grad(x, gradient, N, aux);
    pp0 = vec_dot_prod(gradient, p, N);

    double pi;              /* phi(i)       */
    double ppi;             /* phi'(i)      */
    double pi_1  = p0;      /* phi(i-1)     */
    //double ppi_1 = pp0;     /* phi'(i-1)    */
    double b_i   = beta1;   /* beta_i       */
    double b_i_1 = beta0;   /* beta_{i-1}   */

    int i, j;
    for (i=0;i<steplength_max;i++)
    {
        printf("Starting bracketing, iteration %d\n",i);
        printf("Bracket bounds: %E, %E\n",b_i_1, b_i);
        fflush(stdout);

        /* calculate objective and derivative values at b_i */
        for (j=0;j<N;j++)
        {
            x_update[j] = x[j] + b_i * p[j];
        }
        pi = misfit(x_update, N, aux);

        if ((pi > p0 + WOLFE_C1 * b_i * pp0) || (pi >= pi_1 && i>0))
        {
            b_star = ls_zoom(x, p, N, opt_par, b_i_1, b_i, aux);
            if (b_star < 0)
            {
                /* ls_zoom failed - try bigger steps */
                b_i = 2 * alpha0;
                b_i_1 = alpha0;
                continue;
            }
            success = 1;
            break;
        }

        grad(x_update, gradient, N, aux);
        ppi = vec_dot_prod(gradient, p, N);
        if (fabs(ppi) <= -WOLFE_C2 * pp0)
        {
            b_star = b_i;
            success = 1;
            break;
        }

        if (ppi >= 0.0)
        {
            b_star = ls_zoom(x, p, N, opt_par, b_i, b_i_1, aux);
            if (b_star < 0)
            {
                /* ls_zoom failed - try bigger steps */
                b_i = 2 * alpha0;
                b_i_1 = alpha0;
                continue;
            }
            success = 1;
            break;
        }

        /* update b_i and b_i_1 */
        b_i_1 = b_i;
        b_i = 2 * b_i;

        if (b_i > beta_max)
        {
            fprintf(stderr,"Error: ls_bracket: max step length reached\n");
            success = 0;
            break;
        }
    }

    free(x_update);
    free(gradient);

    if (!success)
    {
        fprintf(stderr,"Line search bracket failed to return an "
        "appropriate step length\n");
        exit(-90);
    }

    printf("Step length successfully selected. Step length: %E\n",b_star);
    fflush(stdout);
    return b_star;
}

/*****************************************************************************/

int check_angle_restart(double *grad, double *p, int N)
{
    int restart = 0;
    double angle_tol_deg = 95.0;
    double angle_tol_rad = angle_tol_deg * (M_PI / 180.0);

    double grad_norm = vec_norm(grad, N);
    double p_norm = vec_norm(p, N);
    double pg = vec_dot_prod(p, grad, N);
    double ang = pg / (p_norm * grad_norm);
    if(ang > cos(angle_tol_rad))
    {
        restart = 1;
    }

    return restart;
}

/*****************************************************************************/

double back_track_quad(double phi0, double phip0, double phi_a0, double a0)
{
    /* see Nodedal and Wright, p. 58 (equation 3.57 and 3.58) */

    double a1;

    double num = phip0 * a0 * a0;
    double den = phi_a0 - phi0 - phip0 * a0;

    a1 = -0.5 * (num / den);
    return a1;
}

double back_track_cubic(double phi0,
                        double phip0,
                        double phi_a0,
                        double phi_a1,
                        double a0,
                        double a1)
{
    /* see Nocedal and Wright 2006, p. 58 */

    double a, b, a2;

    double a02 = a0 * a0;   /* a0 squared   */
    double a12 = a1 * a1;   /* a1 squared   */
    double a03 = a02 * a0;  /* a0 cubed     */
    double a13 = a12 * a1;  /* a1 cubed     */

    double v1 = phi_a1 - phi0 - phip0 * a1;
    double v2 = phi_a0 - phi0 - phip0 * a0;

    double norm = 1.0 / (a02 * a12 * (a1 - a0));

    a = norm * (a02 * v1 - a12 * v2);
    b = norm * (-a03 * v1 + a13 * v2);

    a2 = (-b + sqrt(b * b - 3 * a * phip0)) / (3 * a);

    return a2;
}

/*****************************************************************************/

double back_track(double *x,
                  double *p,
                  int N,
                  optimization_parameters *opt_par,
                  double alpha,
                  void *aux
                  )
{
    static double a[2]; /* holds a set of step lengths that can be used for bracketing */
    a[0] = 0;
    a[1] = 0;

    double (*misfit)(const double *, const int, void *) = opt_par->objective;
    void (*grad)(const double *, double *, int, void *) = opt_par->grad;

    int i;

    double phi0, phip0;
    double a0, a1, a2;
    double phi_a0, phi_a1, phi_a2;

    double *gradient = opt_allocate(sizeof(*gradient) * N);
    double *x_update = opt_allocate(sizeof(*x_update) * N);

    /* evaluate initial misfit and directional derivative */
    phi0 = misfit(x, N, aux);
    grad(x, gradient, N, aux);
    phip0 = vec_dot_prod(gradient, p, N);

    /* try a0*/
    a0 = alpha;
    a[1] = a0;
    for (i=0;i<N;i++)
    {
        x_update[i] = x[i] + a0 * p[i];
    }
    phi_a0 = misfit(x_update, N, aux);

    /* check Armijo condition */
    if (phi_a0 <= phi0 + WOLFE_C1 * a0 * phip0)
    {
        /* condition satisfied */
        free(x_update);
        free(gradient);
        return a0;
    }

    /* try quadratic approximation */
    a1 = back_track_quad(phi0, phip0, phi_a0, a0);
    a[0] = a1;

    for (i=0;i<N;i++)
    {
        x_update[i] = x[i] + a1 * p[i];
    }
    phi_a1 = misfit(x_update, N, aux);

    /* check Armijo condition */
    if (phi_a1 <= phi0 + WOLFE_C1 * a1 * phip0)
    {
        /* condition satisfied */
        free(x_update);
        free(gradient);
        return a1;
    }

    /* move into cubic approximations */
    int cubic_max = 10; /* <----------May need to adjust max number of iterations */
    int j;

    for (j=0;j<cubic_max;j++)
    {
        a2 = back_track_cubic(phi0, phip0, phi_a0, phi_a1, a0, a1);

        /* make sure a2 is not too close to a1 */
        if (fabs(a2 - a1) / a1 < 0.01) /* <--------- May need to adjust the factor of 0.01 here */
            a2 = 0.5 * a1;

        /* make sure a2 is not too much smaller than a1 */
        if (a2 < (0.01 * a1)) /* <-------- May need to adjust the factor of 0.01 here */
            a2 = 0.5 * a1;

        a[1] = a[0];
        a[0] = a2;

        for (i=0;i<N;i++)
        {
            x_update[i] = x[i] + a2 * p[i];
        }
        phi_a2 = misfit(x_update, N, aux);

        /* check Armijo condition */
        if (phi_a2 <= phi0 + WOLFE_C1 * a2 * phip0)
        {
            /* condition satisfied */
            break;
        }

        /* keep only the two newest values */
        a0 = a1;
        phi_a0 = phi_a1;
        a1 = a2;
        phi_a1 = phi_a2;
    }

    if (j >= cubic_max)
    {
        fprintf(stderr,"Backtrack failed to return an "
                       "appropriate step length\n");
        //exit(1);
    }

    free(x_update);
    free(gradient);
    return a2;
}

/*****************************************************************************/

double max3(double x1, double x2, double x3)
{
    double mm = (x1 > x2) ? x1 : x2;
    mm = (mm > x3) ? mm : x3;
    return mm;

}

void modified_chol(double *A, int N)
{

    /* make temp copy of A */
    double *temp = opt_allocate(sizeof(*temp) * N * N);
    double *c = opt_allocate(sizeof(*c) * N * N);
    double *d = opt_allocate(sizeof(*d) * N);
    memcpy(temp, A, sizeof(*temp) * N *N);

    int i, j, s;

    double sum = 0;

    double beta = 1.0;
    double delta = 1.0e-3;
    double theta = 0.0;

    for (j = 0; j < N; j++)
    {
        for (s=0;s<j;s++)
        {
            sum = d[s] * A[j+s*N] * A[j+s*N];
        }

        c[j+j*N] = temp[j+j*N] - sum;
        A[j+j*N] = 1.0;

        //d[j] = c[j+j*N];
        d[j] = max3(fabs(c[j+j*N]),pow(theta / beta,2.0),delta);
        theta = DBL_MIN;
        for (i = j+1; i<N;i++)
        {
            for (s=0;s<j;s++)
            {
                sum = d[s] * A[i+s*N] * A[j+s*N];
            }

            c[i+j*N] = temp[j+i*N] -sum;
            theta = theta > c[i+j*N] ? theta : c[i+j*N];
            A[i+j*N] = c[i+j*N] / d[j];
        }
    }

    for (j=0;j<N;j++)
    {
        for (i=j;i<N;i++)
        {
            A[i+j*N] = A[i+j*N] * sqrt(d[j]);
        }
    }

    //printf("Diag:\n");
    //for (i=0;i<N;i++)
    //    printf("%f\n",d[i]);

    free(d);
    free(c);
    free(temp);
}

void* opt_allocate(long nbytes)
{
    void *ptr = malloc(nbytes);
    if (!ptr)
    {
        fprintf(stderr, "Error: optimization.c : memory allocation failure\n");
        exit (1);
    }
    return ptr;
}

ext: c

max_stars_repo_name: dt-jackson/FWI

max_issues_repo_name: dt-jackson/FWI

max_forks_repo_name: dt-jackson/FWI


text:
// -*- C++ -*-
//
// michael a.g. aïvázis <michael.aivazis@para-sim.com>
//
// (c) 2013-2020 parasim inc
// (c) 2010-2020 california institute of technology
// all rights reserved
//

// code guard
#if !defined(altar_bayesian_COV_h)
#define altar_bayesian_COV_h

// externals
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

// place everything in the local namespace
namespace altar {
    namespace bayesian {

        // forward declarations
        class COV;
        class CoolingStep;

    } // of namespace bayesian
} // of namespace altar

// declaration
class altar::bayesian::COV
{
    // types
public:
    typedef CoolingStep state_t;

    typedef gsl_rng rng_t;
    typedef gsl_vector vector_t;
    typedef gsl_matrix matrix_t;

    // accessors
public:
    inline auto cov() const;
    inline auto beta() const;

    // interface
public:
    virtual double dbeta_brent(vector_t * llk, double llkMedian, vector_t * w);
    virtual double dbeta_grid(vector_t * llk, double llkMedian, vector_t * w);

    // meta-methods
public:
    virtual ~COV();
    inline COV(rng_t * rng, double tolerance=.001, size_t maxIterations=1000, double target=1.0);

    // data
private:
    const double _betaMin;
    const double _betaMax;

    rng_t * _rng;
    double _tolerance;
    size_t _maxIterations;
    double  _target;
    double _beta, _cov;

    // disallow
private:
    COV(const COV &) = delete;
    const COV & operator=(const COV &) = delete;
};

// get the inline definitions
#define altar_bayesian_COV_icc
#include "COV.icc"
#undef altar_bayesian_COV_icc

# endif
// end of file

max_stars_repo_name: AlTarFramework/altar

max_issues_repo_name: AlTarFramework/altar

max_forks_repo_name: AlTarFramework/altar


text:
#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

/* Level 1 */
int mgsl_blas_caxpy(gsl_complex_float *alpha, gsl_vector_complex_float *x, gsl_vector_complex_float *y)
{
  return gsl_blas_caxpy(*alpha, x, y);
}

int mgsl_blas_zaxpy(gsl_complex *alpha, gsl_vector_complex *x, gsl_vector_complex *y)
{
  return gsl_blas_zaxpy(*alpha, x, y);
}

void mgsl_blas_cscal(gsl_complex_float *alpha, gsl_vector_complex_float *x)
{
  gsl_blas_cscal(*alpha, x);
}

void mgsl_blas_zscal(gsl_complex *alpha, gsl_vector_complex *x)
{
  gsl_blas_zscal(*alpha, x);
}

/* Level 2 */
int mgsl_blas_cgemv(CBLAS_TRANSPOSE_t TransA, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_vector_complex_float *x, gsl_complex_float *beta, gsl_vector_complex_float *y)
{
  return gsl_blas_cgemv(TransA, *alpha, A, x, *beta, y);
}

int mgsl_blas_zgemv(CBLAS_TRANSPOSE_t TransA, gsl_complex *alpha, gsl_matrix_complex *A, gsl_vector_complex *x, gsl_complex *beta, gsl_vector_complex *y)
{
  return gsl_blas_zgemv(TransA, *alpha, A, x, *beta, y);
}

int mgsl_blas_chemv(CBLAS_UPLO_t Uplo, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_vector_complex_float *x, gsl_complex_float *beta, gsl_vector_complex_float *y)
{
  return gsl_blas_chemv(Uplo, *alpha, A, x, *beta, y);
}

int mgsl_blas_zhemv(CBLAS_UPLO_t Uplo, gsl_complex *alpha, gsl_matrix_complex *A, gsl_vector_complex *x, gsl_complex *beta, gsl_vector_complex *y)
{
  return gsl_blas_zhemv(Uplo, *alpha, A, x, *beta, y);
}

int mgsl_blas_cgeru(gsl_complex_float *alpha, gsl_vector_complex_float *x, gsl_vector_complex_float *y, gsl_matrix_complex_float *A)
{
  return gsl_blas_cgeru(*alpha, x, y, A);
}

int mgsl_blas_zgeru(gsl_complex *alpha, gsl_vector_complex *x, gsl_vector_complex *y, gsl_matrix_complex *A)
{
  return gsl_blas_zgeru(*alpha, x, y, A);
}

int mgsl_blas_cgerc(gsl_complex_float *alpha, gsl_vector_complex_float *x, gsl_vector_complex_float *y, gsl_matrix_complex_float *A)
{
  return gsl_blas_cgerc(*alpha, x, y, A);
}

int mgsl_blas_zgerc(gsl_complex *alpha, gsl_vector_complex *x, gsl_vector_complex *y, gsl_matrix_complex *A)
{
  return gsl_blas_zgerc(*alpha, x, y, A);
}

int mgsl_blas_cher2(CBLAS_UPLO_t Uplo, gsl_complex_float *alpha, gsl_vector_complex_float *x, gsl_vector_complex_float *y, gsl_matrix_complex_float *A)
{
  return gsl_blas_cher2(Uplo, *alpha, x, y, A);
}

int mgsl_blas_zher2(CBLAS_UPLO_t Uplo, gsl_complex *alpha, gsl_vector_complex *x, gsl_vector_complex *y, gsl_matrix_complex *A)
{
  return gsl_blas_zher2(Uplo, *alpha, x, y, A);
}

int mgsl_blas_cgemm(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B, gsl_complex_float *beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_cgemm(TransA, TransB, *alpha, A, B, *beta, C);
}

int mgsl_blas_zgemm(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B, gsl_complex *beta, gsl_matrix_complex *C)
{
  return gsl_blas_zgemm(TransA, TransB, *alpha, A, B, *beta, C);
}

int mgsl_blas_csymm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B, gsl_complex_float *beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_csymm(Side, Uplo, *alpha, A, B, *beta, C);
}

int mgsl_blas_zsymm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B, gsl_complex *beta, gsl_matrix_complex *C)
{
  return gsl_blas_zsymm(Side, Uplo, *alpha, A, B, *beta, C);
}

int mgsl_blas_chemm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B, gsl_complex_float *beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_chemm(Side, Uplo, *alpha, A, B, *beta, C);
}

int mgsl_blas_zhemm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B, gsl_complex *beta, gsl_matrix_complex *C)
{
  return gsl_blas_zhemm(Side, Uplo, *alpha, A, B, *beta, C);
}

int mgsl_blas_ctrmm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B)
{
  return gsl_blas_ctrmm(Side, Uplo, TransA, Diag, *alpha, A, B);
}

int mgsl_blas_ztrmm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B)
{
  return gsl_blas_ztrmm(Side, Uplo, TransA, Diag, *alpha, A, B);
}

int mgsl_blas_ctrsm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B)
{
  return gsl_blas_ctrsm(Side, Uplo, TransA, Diag, *alpha, A, B);
}

int mgsl_blas_ztrsm(CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B)
{
  return gsl_blas_ztrsm(Side, Uplo, TransA, Diag, *alpha, A, B);
}

int mgsl_blas_csyrk(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_complex_float *beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_csyrk(Uplo, TransA, *alpha, A, *beta, C);
}

int mgsl_blas_zsyrk(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex *alpha, gsl_matrix_complex *A, gsl_complex *beta, gsl_matrix_complex *C)
{
  return gsl_blas_zsyrk(Uplo, TransA, *alpha, A, *beta, C);
}

int mgsl_blas_csyr2k(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B, gsl_complex_float *beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_csyr2k(Uplo, TransA, *alpha, A, B, *beta, C);
}

int mgsl_blas_zsyr2k(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B, gsl_complex *beta, gsl_matrix_complex *C)
{
  return gsl_blas_zsyr2k(Uplo, TransA, *alpha, A, B, *beta, C);
}

int mgsl_blas_cher2k(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex_float *alpha, gsl_matrix_complex_float *A, gsl_matrix_complex_float *B, float beta, gsl_matrix_complex_float *C)
{
  return gsl_blas_cher2k(Uplo, TransA, *alpha, A, B, beta, C);
}

int mgsl_blas_zher2k(CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA, gsl_complex *alpha, gsl_matrix_complex *A, gsl_matrix_complex *B, double beta, gsl_matrix_complex *C)
{
  return gsl_blas_zher2k(Uplo, TransA, *alpha, A, B, beta, C);
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-BLAS

max_issues_repo_name: frithnanth/raku-Math-Libgsl-BLAS

max_forks_repo_name: frithnanth/raku-Math-Libgsl-BLAS


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

#include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

#include "LSAbundance.h"

static int  verbose  = FALSE;

static int  nLines   = 9;

static char *usage[] = {"LSAbundance - \n",
                        "Required parameters:\n",
			"   -a    abundance data\n",
			"   -in   filename      parameter file  \n",
                        "Optional:\n",
			"   -b    integer       length of sample file to ignore\n",
			"   -s    integer       sampling frequency\n",
			"   -seed long          seed random number generator\n",
			"   -v                  verbose\n"};

void updateExpectations(double* adExpect, int nMax, double dMDash, double dV, double dNu, int nS, t_Data *ptData)
{
  int i = 0;

  for(i = 1; i <= nMax; i++){
    int nA = i;
    double dLog = 0.0, dP = 0.0;

    if(nA < MAX_QUAD){
      dLog = logLikelihoodQuad(nA, dMDash, dV, dNu);
    }
    else{
      dLog = logLikelihoodRampal(nA, dMDash, dV, dNu);
    }

    dP = exp(dLog);

    adExpect[i - 1]+= dP*nS;
  }
}

int main(int argc, char* argv[]){
  t_Params    tParams;
  t_LSParams* atLSParams;
  int         i = 0, nSamples = 0, nMax = 0;
  t_Data      tData;
  double*     adExpect = NULL;

  gsl_set_error_handler_off();

  getCommandLineParams(&tParams, argc, argv);

  /*allocate memory for samples*/
  atLSParams = (t_LSParams *) malloc(MAX_SAMPLES*sizeof(t_LSParams));
  if(!atLSParams)
    goto memoryError;

  /*read in Monte-Carlo samples*/
  readSamples(&tParams, atLSParams, &nSamples);

  readAbundanceData(tParams.szAbundFile, &tData);

  nMax = tData.aanAbund[tData.nNA - 1][0];
  nMax = floor(pow(2.0,ceil(log((double) nMax)/log(2.0)) + 2.0) + 1.0e-7);

  adExpect = (double *) malloc(sizeof(double)*nMax);

  for(i = 0; i < nMax; i++){
    adExpect[i] = 0.0;
  }

  for(i = 0; i < nSamples; i++){
    updateExpectations(adExpect, nMax,
		       atLSParams[i].dMDash,
		       atLSParams[i].dV,
		       atLSParams[i].dNu,
		       atLSParams[i].nS,
		       &tData);
  }

  for(i = 1; i <= nMax; i++){
    printf("%d %f\n",i,adExpect[i - 1]/((double) nSamples));
  }

  free(adExpect);
  exit(EXIT_SUCCESS);

 memoryError:
  fprintf(stderr, "Failed to allocate memory in main aborting ...\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int intCompare(const void *pvA, const void *pvB)
{
  int* pnA = (int *) pvA, *pnB = (int *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

int doubleCompare(const void *pvA, const void *pvB)
{
  double* pnA = (double *) pvA, *pnB = (double *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

 /*get abundance file name*/
  ptParams->szAbundFile  = extractParameter(argc,argv, ABUND_FILE,ALWAYS);
  if(ptParams->szAbundFile == NULL)
    goto error;

  /*get long seed*/
  szTemp  = extractParameter(argc,argv,SEED,OPTION);
  if(szTemp != NULL){
    ptParams->lSeed = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->lSeed = 0;
  }

  /*get burn*/
  szTemp  = extractParameter(argc, argv, BURN, OPTION);
  if(szTemp != NULL){
    ptParams->nBurn = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nBurn = DEF_BURN;
  }

 /*get long seed*/
  szTemp  = extractParameter(argc, argv, SAMPLE, OPTION);
  if(szTemp != NULL){
    ptParams->nSample = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nSample = DEF_SAMPLE;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readSamples(t_Params *ptParams, t_LSParams *atLSParams, int *pnSamples)
{
  int  nSamples = 0;
  char *szInputFile = ptParams->szInputFile;
  char szLine[MAX_LINE_LENGTH];
  FILE *ifp = NULL;

  ifp = fopen(szInputFile, "r");

  if(ifp){
    while(fgets(szLine, MAX_LINE_LENGTH, ifp)){
      char *szTok = NULL, *szBrk = NULL, *pcError = NULL;
      int  nTime  = 0;

      /*remove trailing new line*/
      szBrk = strpbrk(szLine, "\n"); (*szBrk) = '\0';

      szTok = strtok(szLine, DELIM);

      nTime = strtol(szTok, &pcError, 10);
      if(*pcError != '\0') goto fileFormatError;

      if(nTime > ptParams->nBurn && nTime % ptParams->nSample == 0){

	szTok = strtok(NULL,DELIM);
	atLSParams[nSamples].dMDash = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atLSParams[nSamples].dV = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atLSParams[nSamples].dNu = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atLSParams[nSamples].nS = strtol(szTok, &pcError, 10);
	if(*pcError != '\0') goto fileFormatError;

	nSamples++;
      }
    }

    fclose(ifp);
  }
  else{
    fprintf(stderr, "Failed to open file %s aborting\n", szInputFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  (*pnSamples) = nSamples;
  return;

 fileFormatError:
  fprintf(stderr, "Incorrectly formatted input file aborting\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

double f1(double x, void *pvParams)
{
  t_LSParams *ptLSParams = (t_LSParams *) pvParams;
  double dMDash = ptLSParams->dMDash, dV = ptLSParams->dV, dNu = ptLSParams->dNu;
  int n = ptLSParams->n;
  double t = ((x - dMDash)*(x - dMDash))/dV;
  double dExp  = x*((double) n) - exp(x);
  double dF  = pow(1.0 + t/dNu, -0.5*(dNu + 1.0));

  return exp(dExp)*dF;
}

double f1Log(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV;
  int n = ptLNParams->n;
  double dTemp = (x - dMDash);
  double dExp  = x*((double) n) - exp(x) - 0.5*((dTemp*dTemp)/dV);
  double dRet  = exp(dExp);

  return dRet;
}

double derivExponent(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
  double dTemp = (x - dMDash)/dV, dRet = 0.0;

  dRet = ((double) n) - exp(x) - dTemp;

  return dRet;
}

double logStirlingsGamma(double dZ)
{
  return 0.5*log(2.0*M_PI) + (dZ - 0.5)*log(dZ) - dZ;
}

double logLikelihoodQuad(int n, double dMDash, double dV, double dNu)
{
  gsl_integration_workspace *ptGSLWS =
    gsl_integration_workspace_alloc(1000);
  double dLogFac1   = 0.0, dLogFacN  = 0.0;
  double dN = (double) n, dResult = 0.0, dError = 0.0, dPrecision = 0.0;
  gsl_function tGSLF;
  t_LSParams tLSParams;
  double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;

  tLSParams.n = n; tLSParams.dMDash = dMDash; tLSParams.dV = dV; tLSParams.dNu = dNu;

  tGSLF.function = &f1;
  tGSLF.params   = (void *) &tLSParams;

  if(dNu < MAX_MU_GAMMA){
    dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
  }
  else{
    dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
  }

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else if(n < 100){
    dLogFacN = gsl_sf_lngamma(dN + 1.0);
  }
  else{
    dLogFacN = logStirlingsGamma(dN + 1.0);
  }

  dA = -100.0; dB = 100.0;

  if(n < 10){
    dPrecision = HI_PRECISION;
  }
  else{
    dPrecision = LO_PRECISION;
  }

  gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);

  //printf("%f %f\n", dResult, dError);

  gsl_integration_workspace_free(ptGSLWS);

  return log(dResult) - dLogFacN + dLogFac1 - 0.5*log(dV);
}

int solveF(double x_lo, double x_hi, double (*f)(double, void*),
	   void* params, double tol, double *xsolve)
{
  int status, iter = 0, max_iter = 100;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double r = 0;
  gsl_function F;
  t_LNParams *ptLNParams = (t_LNParams *) params;
  F.function = f;
  F.params = params;

  //printf("%f %f %d %f %f\n",ptLNParams->dMDash, ptLNParams->dV, ptLNParams->n, x_lo, x_hi);
  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);

  do{
    iter++;
    status = gsl_root_fsolver_iterate (s);
    r = gsl_root_fsolver_root (s);
    x_lo = gsl_root_fsolver_x_lower (s);
    x_hi = gsl_root_fsolver_x_upper (s);

    status = gsl_root_test_interval (x_lo, x_hi, 0, tol);
  }
  while (status == GSL_CONTINUE && iter < max_iter);

  (*xsolve) = gsl_root_fsolver_root (s);
  gsl_root_fsolver_free (s);

  return status;
}

double logLikelihoodLNQuad(int n, double dMDash, double dV)
{
  gsl_integration_workspace *ptGSLWS =
    gsl_integration_workspace_alloc(1000);
  double dLogFac1   = 0.0, dLogFacN  = 0.0;
  double dResult = 0.0, dError = 0.0, dPrecision = 0.0;
  gsl_function tGSLF;
  t_LNParams tLNParams;
  double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;

  tLNParams.n = n; tLNParams.dMDash = dMDash; tLNParams.dV = dV;

  tGSLF.function = &f1Log;
  tGSLF.params   = (void *) &tLNParams;

  dLogFac1 = log(2.0*M_PI*dV);

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else{
    dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
  }

  if(dEst > dV){
    double dMax = 0.0;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    if(fabs(dUpper) > 1.0e-7){
      solveF(0.0, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);
    }

    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }
  else{
    double dMax = 0.0;
    double dLower = dEst - dV;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    if(fabs(dUpper - dLower) > 1.0e-7){
      solveF(dLower, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);
    }
    else{
      dMax = 0.5*(dLower + dUpper);
    }
    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }

  if(n < 10){
    dPrecision = HI_PRECISION;
  }
  else{
    dPrecision = LO_PRECISION;
  }

  gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);

  gsl_integration_workspace_free(ptGSLWS);

  return log(dResult) - dLogFacN -0.5*dLogFac1;
}

double logLikelihoodLNRampal(int n, double dMDash, double dV)
{
  double dN = (double) n;
  double dLogLik = 0.0, dTemp = gsl_pow_int(log(dN) - dMDash,2), dTemp3 = gsl_pow_int(log(dN) - dMDash,3);

  dLogLik = -0.5*log(2.0*M_PI*dV) - log(dN) - (dTemp/(2.0*dV));

  dLogLik += log(1.0 + 1.0/(2.0*dN*dV)*(dTemp/dV + log(dN) - dMDash - 1.0)
		 + 1.0/(6.0*dN*dN*dV*dV*dV)*(3.0*dV*dV - (3.0*dV - 2.0*dV*dV)*(dMDash - log(dN))
		 - 3.0*dV*dTemp + dTemp3));

  return dLogLik;
}

double logLikelihoodRampal(int n, double dMDash, double dV, double dNu)
{
  double dGamma = 0.5*(dNu + 1.0), dN = (double) n, dRN = 1.0/dN, dRSV = 1.0/(sqrt(dV)*sqrt(dNu));
  double dZ = (log(dN) - dMDash)*dRSV;
  double dDZDX = dRN*dRSV, dDZDX2 = -dRN*dRN*dRSV;
  double dF = (1.0 + dZ*dZ);
  double dA = 0.0, dB = 0.0, dTemp = 0.0;
  double dLogFac1 = 0.0;

  if(dNu < MAX_MU_GAMMA){
    dLogFac1 = gsl_sf_lngamma(0.5*(dNu + 1.0)) - gsl_sf_lngamma(0.5*dNu) - 0.5*log(M_PI*dNu);
  }
  else{
    dLogFac1 = 0.5*dNu*(log(0.5*(dNu + 1.0)) - log(0.5*dNu)) -0.5*log(2.0*M_PI) - 0.5;
  }

  dA = 4.0*dZ*dZ*dDZDX*dDZDX*dGamma*(dGamma + 1.0);
  dA /= dF*dF;

  dB = -2.0*dGamma*(dDZDX*dDZDX + dZ*dDZDX2);
  dB /= dF;

  dTemp = dRN + dA + dB;

  return -dGamma*log(dF) + log(dTemp) + dLogFac1 - 0.5*log(dV);
}

void readAbundanceData(const char *szFile, t_Data *ptData)
{
  int **aanAbund = NULL;
  int  i = 0, nNA = 0, nA = 0, nC = 0;
  int  nL = 0, nJ = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    szTok = strtok(szLine, DELIM2);

    nNA = strtol(szTok,&pcError,10);
    if(*pcError != '\0'){
      goto formatError;
    }

    aanAbund = (int **) malloc(nNA*sizeof(int*));

    for(i = 0; i < nNA; i++){
      aanAbund[i] = (int *) malloc(sizeof(int)*2);

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM2);

      nA = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      szTok = strtok(NULL, DELIM2);

      nC = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      nL += nC;
      nJ += nC*nA;

      aanAbund[i][0]  = nA;
      aanAbund[i][1]  = nC;
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  ptData->nJ          = nJ;
  ptData->nL          = nL;
  ptData->aanAbund    = aanAbund;
  ptData->nNA         = nNA;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


text:
/*
* -----------------------------------------------------------------
*  ODE Solver Library --- ode_lib.c
*  Version: 1.6180
*  Date: Jan 4, 2011
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  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.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the implementation file of a library
*  with ODE solution tools.
* -----------------------------------------------------------------
*/

#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <nvector/nvector_serial.h>
#include <cvode/cvode_dense.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>

#include "../include/ode_lib.h"

/*
*------------------------------------------------------------
*   uround
*
*   This function computes the machine unit roundoff
*   defined as the smallest u such that 1.0 + u > 1.0
*
*   Output:
*   uround - machine unit roundoff
*
*   last update: May 20, 2009
*------------------------------------------------------------
*/

double uround(void)
{
    double u = 1.0;
    double comp = 0.0;

    while(comp != 1.0)
    {
	u *= 0.5;
	comp = 1.0 + u;
    }

    return 2.0*u;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   wnorm
*
*   This function computes v1 and v2 vector weight norm
*   defined as
*   wnorm = sqrt( (1/n)*(sum (v1[i]*v2[i])^2) )
*
*   Input:
*   n     - vectors dimension
*   v     - vector 1
*   v     - vector 2
*
*   Output:
*   wnorm - weight norm
*
*   last update: May 20, 2009
*------------------------------------------------------------
*/

double wnorm(int n,double *v1,double *v2)
{
    int i;
    double wnorm = 0.0;

    for ( i = 0; i < n; i++ )
	wnorm += v1[i]*v1[i]*v2[i]*v2[i];

    return sqrt( wnorm / (double)n );
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   ewset
*
*   This function creates a error weight vector
*   defined as ewv[i] = rtol*abs(v[i]) + atol
*
*   Input:
*   n    - vector dimension
*   v    - vector
*   atol - absolute tolerance
*   rtol - relative tolerance
*
*   Output:
*   ewv  - error weight vector
*
*   last update: May 22, 2009
*------------------------------------------------------------
*/

void ewtset(int n,double *v,double atol,double rtol,double *ewt)
{
    int i;

    for ( i = 0; i < n; i++ )
        ewt[i] =  rtol*fabs(v[i]) + atol;

    return;
}
/*------------------------------------------------------------*/

/*------------------------------------------------------------
*   jacobian
*
*   This function computes the jacobian matrix of ydot = f(y,t)
*   defiend as  Jij = dfi/dyj
*
*   Input:
*   n      - vector dimension
*   f_data - pointer to external data
*   Fy     - reaction mapping
*   y      - composition vector
*   t      - time step
*   atol   - absolute tolerance
*   rtol   - relative tolerance
*
*   Output:
*   J       - jacobian matrix
*   success or error
*
*   last update: Oct 7, 2009
------------------------------------------------------------*/

int jacobian(CVRhsFn f,void *f_data,gsl_vector *Fy,gsl_vector *y,
                    double t,double atol,double rtol,gsl_matrix *J)
{
    unsigned int i, j, N1, N2;
    double roundoff, min_inc_mult;
    double fnorm, minInc, inc, inc_inv, yjsaved, srur;
    double *yd_data  = NULL;
    double *Fyd_data = NULL;
    N_Vector yd      = NULL;
    N_Vector Fyd     = NULL;
    gsl_vector *ewt  = NULL;

    min_inc_mult = 1.0e3;

    /* vectors dimensions */
    N1 = y->size;
    N2 = Fy->size;

    /* error weight vector */
    ewt = gsl_vector_calloc(N1);

    /* computing the unit roundoff */
    roundoff = uround();

    /* memory allocation and startup of yd and Fyd */
    yd  = N_VNew_Serial(N1);
    Fyd = N_VNew_Serial(N2);
    if ( yd == NULL || Fyd == NULL )
        return GSL_ENOMEM;

    /* obtaining yd and Fyd compoments  */
    yd_data  = NV_DATA_S(yd);
    Fyd_data = NV_DATA_S(Fyd);

    /* setting yd equal y vector */
    for ( i = 0; i < N1; i++ )
        yd_data[i] = y->data[i];

    /* setting Fyd equal the null vector */
    for ( i = 0; i < N2; i++ )
        Fyd_data[i] = 0.0;

    /* defing error weight vector */
    ewtset(N1,y->data,atol,rtol,ewt->data);

    /* computing weight norm */
    fnorm = wnorm(N2,Fy->data,ewt->data);

    /* square root of the machine roundoff */
    srur = sqrt(roundoff);

    /* computing disturbance parameter */
    minInc = (fnorm != 0.0) ?
	    (min_inc_mult*fabs(t)*roundoff*((double)N1)*fnorm) : 1.0;

    for ( j = 0; j < N1; j++ )
    {
        /* saving y[j] value */
        yjsaved = yd_data[j];

        /* disturbance */
        inc = GSL_MAX(srur*fabs(yjsaved),minInc/ewt->data[j]);

    	/* disturbing y */
        yd_data[j] += inc;

        /* computing Fy disturbed */
        f(t,yd,Fyd,f_data);

	/* restoring yd[j] original value */
        yd_data[j] = yjsaved;

	/* computing the step */
        inc_inv = 1.0/inc;

        /* computing jacobian matrix column j */
    	for ( i = 0; i < N2; i++ )
            J->data[i*J->tda+j] = inc_inv*(Fyd_data[i] - Fy->data[i]);
    }

    /* releasing allocated memory */
    N_VDestroy_Serial(yd);
    N_VDestroy_Serial(Fyd);
    gsl_vector_free(ewt);
    yd       = NULL;
    Fyd      = NULL;
    yd_data  = NULL;
    Fyd_data = NULL;
    ewt      = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*------------------------------------------------------------
*   gradient
*
*   This function computes the gradient matrix of R(phi)
*   defiend as: Aij = d R_i/d phi_j
*
*   Input:
*   f       - right hand side function
*   f_data  - pointer to external data
*   t0      - initial time
*   delta_t - time step
*   atol    - absolute tolerance
*   rtol    - relative tolerance
*   phi     - composition vector
*   Rphi    - reac tion mapping
*
*   Output:
*   A       - gradient matrix
*   success or error
*
*   last update: Feb 19, 2010
------------------------------------------------------------*/

int gradient(CVRhsFn f,void *f_data,void *cvode_mem,double t0,double delta_t,
    double atol,double rtol,gsl_vector *phi,gsl_vector *Rphi,gsl_matrix *A)
{
    unsigned int i, j, flag;
    double inc, inc_inv, phijsaved, srur;
    gsl_vector *Rphid = NULL;

    /* memory allocation for Rphid */
    Rphid = gsl_vector_calloc(Rphi->size);

    /* square root of 1.0e6 times the machine roundoff */
    srur = sqrt(1.0e6*DBL_EPSILON);

    for ( j = 0; j < phi->size; j++ )
    {
        /* saving phi_j value */
        phijsaved = phi->data[j];

        /* disturbance */
        inc = phijsaved*srur + srur;

    	/* disturbing phi_j */
        phi->data[j] += inc;

        /* computing R(phi) disturbed */
	flag = odesolver_reinit(f,f_data,t0,atol,rtol,phi,cvode_mem);
        if ( flag != GSL_SUCCESS )
            return flag;

	flag = odesolver(cvode_mem,delta_t,Rphid);
        if ( flag != GSL_SUCCESS )
            return flag;

	/* restoring phi_j original value */
        phi->data[j] = phijsaved;

	/* computing the step */
        inc_inv = 1.0/inc;

        /* computing gradient matrix column j */
    	for ( i = 0; i < Rphi->size; i++ )
            A->data[i*A->tda+j] = inc_inv*(Rphid->data[i] - Rphi->data[i]);
    }

    /* releasing allocated memory */
    gsl_vector_free(Rphid);
    Rphid = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   linear_approx
*
*   This function computes the linear approximation
*   for a vector function F(x) near the point x0.
*
*   F(x) = F(x0) + DF(x0)*(x-x0) + O(|x-x0|^2)
*
*   Input:
*   x    - vector x
*   x0   - vector x0
*   Fx0  - function at x0
*   DFx0 - jacobian matrix at x0
*
*   Output:
*   Fx   - linear approximation for F(x)
*
*   last update: Jun 9, 2009
*------------------------------------------------------------
*/

void linear_approx(gsl_vector *x,gsl_vector *x0,
                gsl_vector *Fx0,gsl_matrix *DFx0,gsl_vector *Fx)
{
    gsl_vector *aux = NULL;

    /* memory allocation */
    aux = gsl_vector_alloc(x0->size);

    /* aux := x */
    gsl_vector_memcpy(aux,x);

    /* aux := x - x0 */
    gsl_vector_sub(aux,x0);

    /* Fx := F(x0) */
    gsl_vector_memcpy(Fx,Fx0);

    /* Fx := F(x0) + DF(x0)*(x-x0) */
    gsl_blas_dgemv(CblasNoTrans,1.0,DFx0,aux,1.0,Fx);

    /* releasing allocated memory */
    gsl_vector_free(aux);
    aux = NULL;

    return;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   odesolver_init
*
*   This function initiates the ODE solver workspace.
*
*   Input:
*   f         - right hand side function
*   f_data    - right hand side function data
*   t0        - initial time
*   mxsteps   - max # of solver iterations
*   atol      - absolute tolerance
*   rtol      - relative tolerance
*   x         - initial condition
*   cvode_mem - ODE solver workspace
*
*   Output:
*   success or error
*
*   last update: Jan 4, 2011
*------------------------------------------------------------
*/

int odesolver_init(CVRhsFn f,void *f_data,double t0,int mxsteps,
		   double atol,double rtol,gsl_vector *x,void *cvode_mem)
{
    int maxnef  = 100;
    int flag    = CV_SUCCESS;
    N_Vector y0 = NULL;

    /* memory allocation and startup of y0 */
    y0 = N_VMake_Serial(x->size,x->data);
    if ( y0 == NULL )
	return GSL_ENOMEM;

    /* memory allocation for CVODE */
    flag = CVodeMalloc(cvode_mem,f,t0,y0,CV_SS,rtol,&atol);
    if( flag != CV_SUCCESS )
    	return GSL_ENOMEM;

    /* setting user data for right hand side function */
    flag = CVodeSetFdata(cvode_mem,f_data);
    if( flag != CV_SUCCESS )
        return flag;

    /* setting max # of steps to be taken by the solver */
    flag = CVodeSetMaxNumSteps(cvode_mem,mxsteps);
    if( flag != CV_SUCCESS )
        return flag;

    /* setting max # of error test failures permitted */
    flag = CVodeSetMaxErrTestFails(cvode_mem,maxnef);
    if( flag != CV_SUCCESS )
        return flag;

    /* setting dense linear solver */
    flag = CVDense(cvode_mem,x->size);
    if( flag != CV_SUCCESS )
        return flag;

    /* releasing allocated memory */
    N_VDestroy_Serial(y0);
    y0 = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   odesolver_reinit
*
*   This function reinitiates the ODE solver workspace.
*
*   Input:
*   f         - right hand side function
*   f_data    - right hand side function data
*   t0        - initial time
*   atol      - absolute tolerance
*   rtol      - relative tolerance
*   x         - initial condition
*   cvode_mem - ODE solver workspace
*
*   Output:
*   success or error
*
*   last update: Oct 31, 2009
*------------------------------------------------------------
*/

int odesolver_reinit(CVRhsFn f,void *f_data,double t0,
		   double atol,double rtol,gsl_vector *x,void *cvode_mem)
{
    int flag    = CV_SUCCESS;
    N_Vector y0 = NULL;

    /* memory allocation and startup of y0 */
    y0 = N_VMake_Serial(x->size,x->data);
    if ( y0 == NULL )
	return GSL_ENOMEM;

    /* restarting CVODE workspace */
    flag = CVodeReInit(cvode_mem,f,t0,y0,CV_SS,rtol,&atol);
    if( flag != CV_SUCCESS )
    	return flag;

    /* setting user data for right hand side function */
    flag = CVodeSetFdata(cvode_mem,f_data);
    if( flag != CV_SUCCESS )
        return flag;

    /* releasing allocated memory */
    N_VDestroy_Serial(y0);
    y0 = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   odesolver
*
*   This function performs the direct integration of the
*   governing equations.
*
*   Input:
*   cvode_mem - ODE solver workspace
*   tf   - final time
*   Fx   - solution vector
*
*   Output:
*   success or error
*
*   last update: Oct 31, 2009
*------------------------------------------------------------
*/

int odesolver(void *cvode_mem,double tf,gsl_vector *Fx)
{
    double t;
    int flag   = CV_SUCCESS;
    N_Vector y = NULL;

    /* memory allocation and startup of y */
    y = N_VMake_Serial(Fx->size,Fx->data);
    if ( y == NULL )
	return GSL_ENOMEM;

    /* calling CVode solver */
    flag = CVode(cvode_mem,tf,y,&t,CV_NORMAL);
    if( flag != CV_SUCCESS )
        return flag;

    /* releasing allocated memory */
    N_VDestroy_Serial(y);
    y = NULL;

    return GSL_SUCCESS;
}
/*------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib