/* This program generates a random series of pitches, using key-profiles, a pitch proximity profile,
and a range profile. Using param_set = 2, one can generate melodies using the generative model
assumed in ch. 4 of _Music and Probability_. 

To reproduce melodies in book:
   
In all cases, mean pitch = 68, key = C. 

param_set = 0 (flat key profile and flat prox profiles; range var = 10.6; seed=1) generates
figure 4.5A.

param_set = 1 (normal prox and range profiles; range var = 29.0; flat key profile; seed=1)
- not used in book.

param_set = 2 (normal key, prox, and range profiles; range var = 29.0; seed=4) - generates
figure 4.5B. This is essentially the 

*/

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

struct {
    int ontime;
    int offtime;
    int pitch;
    int pc;
} note[20000];

int numnotes;

int param_set;

/* Essen profiles */
double major_key_profile[] =  {0.184,  0.001,  0.155,  0.003,  0.191,  0.109,  0.005,  0.214,  0.001,  0.078,  0.004,  0.055};
double minor_key_profile[] =  {0.192,  0.005,  0.149,  0.179,  0.002,  0.144,  0.002,  0.201,  0.038,  0.012,  0.053,  0.022};

/* Theoretical profiles */
/*
double major_key_profile[] = {0.17, 0.01, 0.12, 0.01, 0.15, 0.12, 0.01, 0.15, 0.01, 0.12, 0.01, 0.12};
double minor_key_profile[] = {0.17, 0.01, 0.12, 0.15, 0.01, 0.12, 0.01, 0.15, 0.12, 0.01, 0.01, 0.12};
*/

double mean_profile[100];
double proximity_profile[100];
double range_profile[100];
double adjusted_profile[100];
    
int abs(int x) {
    if(x>=0) return x;
    else return 0-x;
}

int verbosity;

main(argc, argv) 
     int argc;
     char *argv[];
{

    int i, m, seed, mean;
    double total_prox_mass=0.0;
    double var, x;

    param_set = 1;

    /* choose the seed from the time... */
    seed = time(NULL) % 10000;
    /* ...or set the seed here */
    if(param_set == 0 || param_set == 1) seed = 1;
    if(param_set == 2) seed = 4;

    srand(seed);

    verbosity=1;

    /* Set melody length here */
    numnotes = 30;

    x = (double)(rand() % 10000) / 10000.0;
    mean = choose_mean(x);

    if(verbosity>0) printf("Mean = %d\n", mean);

    /* Make proximity profile */
    var = 7.2; /* derived from Essen data */
    for(i=0; i<100; i++) {
	proximity_profile[i] = (exp( -pow( ((double)(i)-50.0), 2.0) / (2.0 * var))) / (2.51 * sqrt(var));
	/* f(x), where x is an interval of 0: f(0) = 1.0; f(1) = e ^ -(.04); f(10) = e ^ -4 = .02 */

	/* To generate flat proximity profile: */
	if(param_set==0) proximity_profile[i] = 1.0;

    }

    /* Make range profile */
    if(param_set==0) var = 10.6; /* crude estimate from Essen data */
    else if(param_set==1 || param_set==2) var = 29.0; /* better estimate from Essen data */
    for(i=0; i<100; i++) {
	range_profile[i] = (exp( -pow( ((double)(i)-50.0), 2.0) / (2.0 * var))) / (2.51 * sqrt(var));
	/* To generate flat range profile: */
	//range_profile[i] = 1.0;
    }

    /* To flatten key profiles: */
    if(param_set==0 || param_set==1) {
	for(i=0; i<12; i++) {
	    major_key_profile[i] = 1.0;
	    minor_key_profile[i] = 1.0;
	}
    }

    if(verbosity>1) {
	printf("Proximity profile: \n");
	/* Print out raw proximity profile values */
	//for(i=0; i<100; i++) printf("%6.3f ", proximity_profile[i]);
	//printf("\n"); 

	/* ...Or print out the profile values normalized to sum to 1 */
	for(i=0; i<100; i++) total_prox_mass += proximity_profile[i];
	for(i=38; i<63; i++) printf("%6.3f ", proximity_profile[i] / total_prox_mass);
	printf("\n");
    }

    if(verbosity>0) {
	printf("Sample adjusted profile assuming key of C major, previous note of 60, mean of %d, over range 48-72:\n", mean);
	make_adjusted_profile(60, 0, mean);
	for(i=55; i<=79; i++) printf("%6.3f ", adjusted_profile[i]);
	printf("\n");
    }

    generate_melody(mean);

}

int choose_mean(double randx) {

    double var, mass;
    int i, mean;
    double column[100];

    /* Make mean profile and choose the mean randomly... */

    var = 13.184; /* derived from Essen data */
    mass = 0.0;
    for(i=0; i<100; i++) {
	mean_profile[i] = (exp( -pow( ((double)(i)-50.0), 2.0) / (2.0 * var))) / (2.51 * sqrt(var));
	mass += mean_profile[i];
    }
    for(i=0; i<100; i++) mean_profile[i] = mean_profile[i] / mass;

    column[0] = mean_profile[0];
    for(i=1; i<100; i++) {
	column[i] = column[i-1] + mean_profile[i];
    }
    if(randx < column[0]) mean=18; /* 0+68-50 = 18 */
    else for(i=1; i<100; i++) {
	if(randx > column[i-1] && randx < column[i]) {
	    mean=i+18;
	    break;
	}
    }

    /* ...or just set the mean here */
    mean = 68;
    return mean;
}

make_adjusted_profile(int prev_p, int k, int mean) {

    /* Make a profile for the next note, using the previous pitch, mean pitch, and key. If it's t he first note
       (prev_p = 0), just use the range and key profiles. */

    int i, interval, range_position;
    double raw_profile[100];
    double mass=0.0;
    int top, bottom;

    for(i=0; i<100; i++) {
	interval = i-prev_p;
	if(interval < -50) interval = -50;
	if(interval > 49) interval = 49;
	range_position = i-mean;
	if(range_position < -50) range_position = -50;
	if(range_position > 49) range_position = 49;

	if(prev_p==-1) {
	    if(k<12) raw_profile[i] = major_key_profile[((i+24)-k)%12] * range_profile[range_position+50];
	    else raw_profile[i] = minor_key_profile[((i+24)-k)%12] * range_profile[range_position+50];
	}

	else {
	    if(k<12) raw_profile[i] = major_key_profile[((i+24)-k)%12] * proximity_profile[interval+50]
			 * range_profile[range_position+50];
	    else raw_profile[i] = minor_key_profile[((i+24)-k)%12] * proximity_profile[interval+50]
		     * range_profile[range_position+50];
	}

	/* To exclude everything outside middle pitch range */
	//if(i < bottom || i > top) raw_profile[i] = 0.0;

	mass += raw_profile[i];
    }

    /* Now normalize the profile values to sum to 1 */
    for(i=0; i<100; i++) {
	adjusted_profile[i] = raw_profile[i] / mass;
    }
}

make_dumb_adjusted_profile(int p, int k) {

    int i, interval;
    double raw_profile[100];
    double mass=0.0;

    for(i=0; i<100; i++) {
	if(i>47 && i<85) raw_profile[i] = major_key_profile[((i+24)-k)%12];
	else raw_profile[i] = 0.0;
	mass += raw_profile[i];
    }

    /* Now normalize the profile values to sum to 1 */
    for(i=0; i<100; i++) {
	adjusted_profile[i] = raw_profile[i] / mass;
    }

    if(verbosity>1) {
	printf("Adjusted profile\n");
	for(i=0; i<100; i++) printf("%6.3f ", adjusted_profile[i]);
	printf("\n"); 
    }
}

print_pitchname(int p) {

    int pc;
    pc = p%12;
    switch(pc) {
    case 0: printf("C"); break;
    case 1: printf("C#"); break;
    case 2: printf("D"); break;
    case 3: printf("Eb"); break;
    case 4: printf("E"); break;
    case 5: printf("F"); break;
    case 6: printf("F#"); break;
    case 7: printf("G"); break;
    case 8: printf("Ab"); break;
    case 9: printf("A"); break;
    case 10: printf("Bb"); break;
    case 11: printf("B"); break;
    }
    printf("%d ", ((p-pc)/12) - 1);
}

generate_melody(int mean) {

    int i, n, key, p, previous_note;
    double x, mass;
    double column[100];

    /* set the key here; 0 = C major */
    key=0;
    previous_note = -1;

    for(n=0; n<numnotes; n++) {

	make_adjusted_profile(previous_note, key, mean);
	if(verbosity>1) {
	    printf("Adjusted profile\n");
	    for(i=0; i<100; i++) printf("%6.3f ", adjusted_profile[i]);
	    printf("\n"); 
	}
	//make_dumb_adjusted_profile(previous_note, key);

	column[0] = adjusted_profile[0];
	for(i=1; i<100; i++) {
	    column[i] = column[i-1] + adjusted_profile[i];
	}
	if(verbosity>1 && previous_note==-1) {
	    printf("Column:\n");
	    for(i=0; i<100; i++) printf("%6.3f\n", column[i]);
	}
	x = (double)(rand() % 10000) / 10000.0;
	if(x < column[0]) p=0;
	else for(i=1; i<100; i++) {
	    if(x > column[i-1] && x < column[i]) {
		p=i;
		break;
	    }
	}
	if(verbosity>0) {
	    printf("%d(", p);
	    print_pitchname(p);
	    printf(") ");
	}

	//if(previous_note != -1) printf("%d ", p-previous_note);

	previous_note = p;
    }
    printf("\n");
}