/*
    SuperCollider real time audio synthesis system
    Copyright (c) 2002 James McCartney. All rights reserved.
    http://www.audiosynth.com

    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
*/

// Gendyn UGens implemented by Nick Collins

#include "SC_PlugIn.h"

static InterfaceTable* ft;

struct Gendy1 : public Unit // Iannis Xenakis/Marie-Helene Serra GENDYN simulation
{
    double mPhase;
    float mFreqMul, mAmp, mNextAmp, mSpeed, mDur;
    int mMemorySize, mIndex;
    float* mMemoryAmp; // could hard code as 12
    float* mMemoryDur;
};

// following Hoffmann paper from CMJ- primary and secondary random walks
struct Gendy2 : public Unit {
    double mPhase;
    float mFreqMul, mAmp, mNextAmp, mSpeed, mDur;
    int mMemorySize, mIndex;
    float* mMemoryAmp;
    float* mMemoryAmpStep;
    float* mMemoryDur;
    float* mMemoryDurStep;
};


// Random walks as Gendy1 but works out all breakpoints per cycle and normalises time intervals to desired frequency
struct Gendy3 : public Unit {
    double mPhase, mNextPhase, mLastPhase;
    float mSpeed, mFreqMul;
    float mAmp, mNextAmp, mInterpMult;
    int mMemorySize, mIndex;
    float* mMemoryAmp;
    float* mMemoryDur;
    double* mPhaseList;
    float* mAmpList;
};

extern "C" {

void Gendy1_next_k(Gendy1* unit, int inNumSamples);
void Gendy1_Ctor(Gendy1* unit);
void Gendy1_Dtor(Gendy1* unit);

void Gendy2_next_k(Gendy2* unit, int inNumSamples);
void Gendy2_Ctor(Gendy2* unit);
void Gendy2_Dtor(Gendy2* unit);

void Gendy3_next_k(Gendy3* unit, int inNumSamples);
void Gendy3_Ctor(Gendy3* unit);
void Gendy3_Dtor(Gendy3* unit);
}


void Gendy1_Ctor(Gendy1* unit) {
    SETCALC(Gendy1_next_k);

    unit->mFreqMul = unit->mRate->mSampleDur;
    unit->mPhase = 1.0; // should immediately decide on new target
    unit->mAmp = 0.0f;
    unit->mNextAmp = 0.0f;
    unit->mSpeed = 100.f;

    unit->mMemorySize = (int)ZIN0(8); // default is 12
    // printf("memsize %d %f", unit->mMemorySize, ZIN0(8));
    if (unit->mMemorySize < 1)
        unit->mMemorySize = 1;
    unit->mIndex = 0;
    unit->mMemoryAmp = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    unit->mMemoryDur = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    ClearUnitIfMemFailed(unit->mMemoryAmp && unit->mMemoryDur);

    RGen& rgen = *unit->mParent->mRGen;

    // initialise to zeroes and separations
    for (int i = 0; i < unit->mMemorySize; ++i) {
        unit->mMemoryAmp[i] = 2 * rgen.frand() - 1.0f;
        unit->mMemoryDur[i] = rgen.frand();
    }

    // compute one sample of output to avoid sending garbage memory downstream to other Ctor functions
    // first sample of the _next output will be the current amplitude (which is 0)
    OUT0(0) = 0.0f;
}

void Gendy1_Dtor(Gendy1* unit) {
    RTFree(unit->mWorld, unit->mMemoryAmp);
    RTFree(unit->mWorld, unit->mMemoryDur);
}


// called once per period so OK to work out constants in here
static float Gendyn_distribution(int which, float a, float f) {
    float temp, c;

    if (a > 1.f)
        a = 1.f; // a must be in range 0 to 1
    if (a < 0.0001f)
        a = 0.0001f; // for safety with some distributions, don't want divide by zero errors

    switch (which) {
    case 0: // LINEAR
            // linear
        break;
    case 1: // CAUCHY
        // X has a*tan((z-0.5)*pi)
        // I went back to first principles of the Cauchy distribution and re-integrated with a
        // normalisation constant

        // choice of 10 here is such that f=0.95 gives about 0.35 for temp, could go with 2 to make it finer
        c = atan(10 * a); // PERHAPS CHANGE TO a=1/a;
        // incorrect- missed out divisor of pi in norm temp= a*tan(c*(2*pi*f - 1));
        temp = (1.f / a) * tan(c * (2.f * f - 1.f)); // Cauchy distribution, C is precalculated

        // printf("cauchy f %f c %f temp %f out %f \n",f,  c, temp, temp/10);

        return temp * 0.1f; //(temp+100)/200;

    case 2: // LOGIST (ic)
        // X has -(log((1-z)/z)+b)/a which is not very usable as is

        c = 0.5f + (0.499f * a); // calculate normalisation constant
        c = log((1.f - c) / c);

        // remap into range of valid inputs to avoid infinities in the log

        // f= ((f-0.5)*0.499*a)+0.5;
        f = ((f - 0.5f) * 0.998f * a) + 0.5f; //[0,1]->[0.001,0.999]; squashed around midpoint 0.5 by a
        // Xenakis calls this the LOGIST map, it's from the range [0,1] to [inf,0] where 0.5->1
        // than take natural log. to avoid infinities in practise I take [0,1] -> [0.001,0.999]->[6.9,-6.9]
        // an interesting property is that 0.5-e is the reciprocal of 0.5+e under (1-f)/f
        // and hence the logs are the negative of each other
        temp = log((1.f - f) / f) / c; // n range [-1,1]
        // X also had two constants in his- I don't bother

        // printf("logist f %f temp %f\n", f, temp);

        return temp; // a*0.5*(temp+1.0);	//to [0,1]

    case 3: // HYPERBCOS
        // X original a*log(tan(z*pi/2)) which is [0,1]->[0,pi/2]->[0,inf]->[-inf,inf]
        // unmanageable in this pure form
        c = tan(1.5692255f * a); // tan(0.999*a*pi*0.5);    	//[0, 636.6] maximum range
        temp = tan(1.5692255f * a * f) / c; //[0,1]->[0,1]
        temp = log(temp * 0.999f + 0.001f) * (-0.1447648f); // multiplier same as /(-6.9077553); //[0,1]->[0,1]

        // printf("hyperbcos f %f c %f temp %f\n", f, c, temp);

        return 2.f * temp - 1.0f;

    case 4: // ARCSINE
        // X original a/2*(1-sin((0.5-z)*pi)) aha almost a better behaved one though [0,1]->[2,0]->[a,0]
        c = sin(1.5707963f * a); // sin(pi*0.5*a);	//a as scaling factor of domain of sine input to use
        temp = sin(pi_f * (f - 0.5f) * a) / c; //[-1,1] which is what I need

        // printf("arcsine f %f c %f temp %f\n", f, c, temp);

        return temp;

    case 5: // EXPON
        // X original -(log(1-z))/a [0,1]-> [1,0]-> [0,-inf]->[0,inf]
        c = log(1.f - (0.999f * a));
        temp = log(1.f - (f * 0.999f * a)) / c;

        // printf("expon f %f c %f temp %f\n", f, c, temp);

        return 2.f * temp - 1.f;

    case 6: // SINUS
        // X original a*sin(smp * 2*pi/44100 * b) ie depends on a second oscillator's value-
        // hmmm, plug this in as a I guess, will automatically accept control rate inputs then!
        return 2.f * a - 1.f;

    default:
        break;
    }

    return 2.f * f - 1.f;
}


void Gendy1_next_k(Gendy1* unit, int inNumSamples) {
    float* out = ZOUT(0);

    // distribution choices for amp and dur and constants of distribution
    int whichamp = (int)ZIN0(0);
    int whichdur = (int)ZIN0(1);
    float aamp = ZIN0(2);
    float adur = ZIN0(3);
    float minfreq = ZIN0(4);
    float maxfreq = ZIN0(5);
    float scaleamp = ZIN0(6);
    float scaledur = ZIN0(7);

    float rate = unit->mDur;

    // phase gives proportion for linear interpolation automatically
    double phase = unit->mPhase;

    float amp = unit->mAmp;
    float nextamp = unit->mNextAmp;

    float speed = unit->mSpeed;

    RGen& rgen = *unit->mParent->mRGen;
    // linear distribution 0.0 to 1.0 using rgen.frand()

    LOOP1(
        inNumSamples, float z;

        if (phase >= 1.0) {
            phase -= 1.0;

            int index = unit->mIndex;
            int num = (int)(ZIN0(9)); //(unit->mMemorySize);(((int)ZIN0(9))%(unit->mMemorySize))+1;
            if ((num > (unit->mMemorySize)) || (num < 1))
                num = unit->mMemorySize;

            // new code for indexing
            index = (index + 1) % num;
            amp = nextamp;

            unit->mIndex = index;

            // Gendy dist gives value [-1,1], then use scaleamp
            // first term was amp before, now must check new memory slot
            nextamp = (unit->mMemoryAmp[index]) + (scaleamp * Gendyn_distribution(whichamp, aamp, rgen.frand()));

            // mirroring for bounds- safe version
            if (nextamp > 1.0f || nextamp < -1.0f) {
                // printf("mirroring nextamp %f ", nextamp);

                // to force mirroring to be sensible
                if (nextamp < 0.0f)
                    nextamp = nextamp + 4.f;

                nextamp = fmod(nextamp, 4.f);
                // printf("fmod  %f ", nextamp);

                if (nextamp > 1.0f && nextamp < 3.f)
                    nextamp = 2.f - nextamp;

                else if (nextamp > 1.0f)
                    nextamp = nextamp - 4.f;

                // printf("mirrorednextamp %f \n", nextamp);
            };

            unit->mMemoryAmp[index] = nextamp;

            // Gendy dist gives value [-1,1]
            rate = (unit->mMemoryDur[index]) + (scaledur * Gendyn_distribution(whichdur, adur, rgen.frand()));

            if (rate > 1.0f || rate < 0.0f) {
                if (rate < 0.0)
                    rate = rate + 2.f;
                rate = fmod(rate, 2.0f);
                rate = 2.f - rate;
            }

            unit->mMemoryDur[index] = rate;

            // printf("nextamp %f rate %f \n", nextamp, rate);

            // define range of speeds (say between 20 and 1000 Hz)
            // can have bounds as fourth and fifth inputs
            speed = (minfreq + ((maxfreq - minfreq) * rate)) * (unit->mFreqMul);

            // if there are 12 control points in memory, that is 12 per cycle
            // the speed is multiplied by 12
            //(I don't store this because updating rates must remain in range [0,1]
            speed *= num;
        }

        // linear interpolation could be changed
        z = ((1.0 - phase) * amp) + (phase * nextamp);

        phase += speed; ZXP(out) = z;);

    unit->mPhase = phase;
    unit->mAmp = amp;
    unit->mNextAmp = nextamp;
    unit->mSpeed = speed;
    unit->mDur = rate;
}


void Gendy2_Ctor(Gendy2* unit) {
    SETCALC(Gendy2_next_k);

    unit->mFreqMul = unit->mRate->mSampleDur;
    unit->mPhase = 1.0; // should immediately decide on new target
    unit->mAmp = 0.0f;
    unit->mNextAmp = 0.0f;
    unit->mSpeed = 100.f;

    unit->mMemorySize = (int)ZIN0(8); // default is 12
    // printf("memsize %d %f", unit->mMemorySize, ZIN0(8));
    if (unit->mMemorySize < 1)
        unit->mMemorySize = 1;
    unit->mIndex = 0;
    unit->mMemoryAmp = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    unit->mMemoryDur = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    unit->mMemoryAmpStep = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    unit->mMemoryDurStep = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    ClearUnitIfMemFailed(unit->mMemoryAmp && unit->mMemoryDur && unit->mMemoryAmpStep && unit->mMemoryDurStep);

    RGen& rgen = *unit->mParent->mRGen;

    // initialise to zeroes and separations
    for (int i = 0; i < unit->mMemorySize; ++i) {
        unit->mMemoryAmp[i] = 2 * rgen.frand() - 1.0f;
        unit->mMemoryDur[i] = rgen.frand();
        unit->mMemoryAmpStep[i] = 2 * rgen.frand() - 1.0f;
        unit->mMemoryDurStep[i] = 2 * rgen.frand() - 1.0f;
    }

    // compute one sample of output to avoid sending garbage memory downstream to other Ctor functions
    // first sample of the _next output will be the current amplitude (which is 0)
    OUT0(0) = 0.0f;
}

void Gendy2_Dtor(Gendy2* unit) {
    RTFree(unit->mWorld, unit->mMemoryAmp);
    RTFree(unit->mWorld, unit->mMemoryDur);
    RTFree(unit->mWorld, unit->mMemoryAmpStep);
    RTFree(unit->mWorld, unit->mMemoryDurStep);
}

static float Gendyn_mirroring(float lower, float upper, float in) {
    // mirroring for bounds- safe version
    if (in > upper || in < lower) {
        float range = (upper - lower);

        if (in < lower)
            in = (2.0f * upper - lower) - in;

        in = fmod(in - upper, 2.0f * range);

        if (in < range)
            in = upper - in;
        else
            in = in - (range);
    }

    return in;
}

void Gendy2_next_k(Gendy2* unit, int inNumSamples) {
    float* out = ZOUT(0);

    // distribution choices for amp and dur and constants of distribution
    int whichamp = (int)ZIN0(0);
    int whichdur = (int)ZIN0(1);
    float aamp = ZIN0(2);
    float adur = ZIN0(3);
    float minfreq = ZIN0(4);
    float maxfreq = ZIN0(5);
    float scaleamp = ZIN0(6);
    float scaledur = ZIN0(7);

    float rate = unit->mDur;

    // phase gives proportion for linear interpolation automatically
    double phase = unit->mPhase;

    float amp = unit->mAmp;
    float nextamp = unit->mNextAmp;

    float speed = unit->mSpeed;

    RGen& rgen = *unit->mParent->mRGen;

    LOOP1(
        inNumSamples, float z;

        if (phase >= 1.0) {
            phase -= 1.0;

            int index = unit->mIndex;
            int num = (int)(ZIN0(9)); //(unit->mMemorySize);(((int)ZIN0(9))%(unit->mMemorySize))+1;
            if ((num > (unit->mMemorySize)) || (num < 1))
                num = unit->mMemorySize;

            // new code for indexing
            index = (index + 1) % num;

            // using last amp value as seed
            // random values made using a lehmer number generator xenakis style
            float a = ZIN0(10);
            float c = ZIN0(11);

            float lehmerxen = fmod(((amp)*a) + c, 1.0f);

            // printf("lehmer %f \n", lehmerxen);

            amp = nextamp;

            unit->mIndex = index;

            // Gendy dist gives value [-1,1], then use scaleamp
            // first term was amp before, now must check new memory slot

            float ampstep = (unit->mMemoryAmpStep[index]) + Gendyn_distribution(whichamp, aamp, fabs(lehmerxen));
            ampstep = Gendyn_mirroring(-1.0f, 1.0f, ampstep);

            unit->mMemoryAmpStep[index] = ampstep;

            nextamp = (unit->mMemoryAmp[index]) + (scaleamp * ampstep);

            nextamp = Gendyn_mirroring(-1.0f, 1.0f, nextamp);

            unit->mMemoryAmp[index] = nextamp;

            float durstep = (unit->mMemoryDurStep[index]) + Gendyn_distribution(whichdur, adur, rgen.frand());
            durstep = Gendyn_mirroring(-1.0f, 1.0f, durstep);

            unit->mMemoryDurStep[index] = durstep;

            rate = (unit->mMemoryDur[index]) + (scaledur * durstep);

            rate = Gendyn_mirroring(0.0f, 1.0f, rate);

            unit->mMemoryDur[index] = rate;

            // printf("nextamp %f rate %f \n", nextamp, rate);

            // define range of speeds (say between 20 and 1000 Hz)
            // can have bounds as fourth and fifth inputs
            speed = (minfreq + ((maxfreq - minfreq) * rate)) * (unit->mFreqMul);

            // if there are 12 control points in memory, that is 12 per cycle
            // the speed is multiplied by 12
            //(I don't store this because updating rates must remain in range [0,1]
            speed *= num;
        }

        // linear interpolation could be changed
        z = ((1.0 - phase) * amp) + (phase * nextamp);

        phase += speed; ZXP(out) = z;);

    unit->mPhase = phase;
    unit->mAmp = amp;
    unit->mNextAmp = nextamp;
    unit->mSpeed = speed;
    unit->mDur = rate;
}

void Gendy3_Ctor(Gendy3* unit) {
    SETCALC(Gendy3_next_k);

    unit->mFreqMul = unit->mRate->mSampleDur;
    unit->mPhase = 1.0; // should immediately decide on new target
    unit->mAmp = 0.0f;
    unit->mNextAmp = 0.0f;
    unit->mNextPhase = 0.0;
    unit->mLastPhase = 0.0;
    unit->mInterpMult = 1.0f;
    unit->mSpeed = 100.f;

    unit->mMemorySize = (int)ZIN0(7);
    if (unit->mMemorySize < 1)
        unit->mMemorySize = 1;
    unit->mIndex = 0;
    unit->mMemoryAmp = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));
    unit->mMemoryDur = (float*)RTAlloc(unit->mWorld, unit->mMemorySize * sizeof(float));

    // one more in amp list for guard (wrap) element
    unit->mAmpList = (float*)RTAlloc(unit->mWorld, (unit->mMemorySize + 1) * sizeof(float));
    unit->mPhaseList = (double*)RTAlloc(unit->mWorld, (unit->mMemorySize + 1) * sizeof(double));
    ClearUnitIfMemFailed(unit->mMemoryAmp && unit->mMemoryDur && unit->mAmpList && unit->mPhaseList);

    RGen& rgen = *unit->mParent->mRGen;

    // initialise to zeroes and separations
    for (int i = 0; i < unit->mMemorySize; ++i) {
        unit->mMemoryAmp[i] = 2 * rgen.frand() - 1.0f;
        unit->mMemoryDur[i] = rgen.frand();
        unit->mAmpList[i] = 2 * rgen.frand() - 1.0f;
        unit->mPhaseList[i] = 1.0; // will be intialised immediately
    }

    unit->mMemoryAmp[0] = 0.0f; // always zeroed first BP

    // compute one sample of output to avoid sending garbage memory downstream to other Ctor functions
    // first sample of the _next output will be the current amplitude (which is 0)
    OUT0(0) = 0.0f;
}

void Gendy3_Dtor(Gendy3* unit) {
    RTFree(unit->mWorld, unit->mMemoryAmp);
    RTFree(unit->mWorld, unit->mMemoryDur);
    RTFree(unit->mWorld, unit->mAmpList);
    RTFree(unit->mWorld, unit->mPhaseList);
}

void Gendy3_next_k(Gendy3* unit, int inNumSamples) {
    float* out = ZOUT(0);

    // distribution choices for amp and dur and constants of distribution
    int whichamp = (int)ZIN0(0);
    int whichdur = (int)ZIN0(1);
    float aamp = ZIN0(2);
    float adur = ZIN0(3);
    float freq = ZIN0(4);
    float scaleamp = ZIN0(5);
    float scaledur = ZIN0(6);

    double phase = unit->mPhase;
    float amp = unit->mAmp;
    float nextamp = unit->mNextAmp;
    float speed = unit->mSpeed;
    int index = unit->mIndex;
    int interpmult = (int)unit->mInterpMult;
    double lastphase = unit->mLastPhase;
    double nextphase = unit->mNextPhase;

    RGen& rgen = *unit->mParent->mRGen;

    float* amplist = unit->mAmpList;
    double* phaselist = unit->mPhaseList;

    LOOP1(
        inNumSamples,
        float z;

        if (phase >= 1.) { // calculate all targets for new period
            phase -= 1.;

            int num = (int)(ZIN0(8));
            if ((num > (unit->mMemorySize)) || (num < 1))
                num = unit->mMemorySize;

            float dursum = 0.0f;

            float* memoryamp = unit->mMemoryAmp;
            float* memorydur = unit->mMemoryDur;

            for (int j = 0; j < num; ++j) {
                if (j > 0) { // first BP always stays at 0
                    float amp = (memoryamp[j]) + (scaleamp * Gendyn_distribution(whichamp, aamp, rgen.frand()));
                    amp = Gendyn_mirroring(-1.0f, 1.0f, amp);
                    memoryamp[j] = amp;
                }

                float dur = (memorydur[j]) + (scaledur * Gendyn_distribution(whichdur, adur, rgen.frand()));
                dur = Gendyn_mirroring(0.01f, 1.0f, dur); // will get normalised in a moment, don't allow zeroes
                memorydur[j] = dur;
                dursum += dur;
            }

            // normalising constant
            dursum = 1.f / dursum;

            int active = 0;

            // phase duration of a sample
            float minphase = unit->mFreqMul;

            speed = freq * minphase;

            // normalise and discard any too short (even first)
            for (int j = 0; j < num; ++j) {
                float dur = memorydur[j];
                dur *= dursum;

                if (dur >= minphase) {
                    amplist[active] = memoryamp[j];
                    phaselist[active] = dur;
                    ++active;
                }
            }

            // add a zero on the end at active
            amplist[active] = 0.0f; // guard element
            phaselist[active] = 2.0; // safety element

            // lastphase=0.0;
            //					nextphase= phaselist[0];
            //					amp=amplist[0];
            //					nextamp=amplist[1];
            //					index=0;
            //					unit->mIndex=index;
            //
            // setup to trigger next block
            nextphase = 0.0;
            nextamp = amplist[0];
            index = -1;
        }

        if (phase >= nextphase) { // are we into a new region?
            // new code for indexing
            ++index; //=index+1; //%num;

            amp = nextamp;

            unit->mIndex = index;

            lastphase = nextphase;
            nextphase = lastphase + phaselist[index];
            nextamp = amplist[index + 1];

            interpmult = (int)(1.0 / (nextphase - lastphase));
        }

        float interp = (phase - lastphase) * interpmult;

        // linear interpolation could be changed
        z = ((1.0f - interp) * amp) + (interp * nextamp);

        phase += speed; ZXP(out) = z;);

    unit->mPhase = phase;
    unit->mSpeed = speed;
    unit->mInterpMult = interpmult;
    unit->mAmp = amp;
    unit->mNextAmp = nextamp;
    unit->mLastPhase = lastphase;
    unit->mNextPhase = nextphase;
}

PluginLoad(Gendyn) {
    ft = inTable;

    DefineDtorUnit(Gendy1);
    DefineDtorUnit(Gendy2);
    DefineDtorUnit(Gendy3);
}