C#，数值计算——分类与推理Phylagglom的计算方法与源程序

1 文本格式

using System;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
public abstract class Phylagglom
{
public int n { get; set; }
public int root { get; set; }
public int fsroot { get; set; }
public double seqmax { get; set; }
public double depmax { get; set; }
public List<Phylagglomnode> t { get; set; } = new List<Phylagglomnode>();
public abstract void premin(double[,] d, int[] nextp);
public abstract double dminfn(double[,] d, int i, int j);
public abstract double dbranchfn(double[,] d, int i, int j);
public abstract double dnewfn(double[,] d, int k, int i, int j, int ni, int nj);
public abstract void drootbranchfn(double[,] d, int i, int j, int ni, int nj, ref double bi, ref double bj);

public Phylagglom(double[,] dist, int fsr = -1)
{
this.n = dist.GetLength(0);
this.fsroot = fsr;
this.t = new List<Phylagglomnode>(2 * n - 1);
}

public void makethetree(double[,] dist)
{
double[,] d = Globals.CopyFrom(dist);
int[] tp = new int[n];
int[] nextp = new int[n];
int[] prevp = new int[n];
int[] tasklist = new int[2 * n + 1];
double[] tmp = new double[n];

int j;
int i = 0;
for (; i < n; i++)
{
nextp[i] = i + 1;
prevp[i] = i - 1;
tp[i] = i;
t[i].ldau = t[i].rdau = -1;
t[i].nel = 1;
}

prevp[0] = nextp[n - 1] = -1;
int ncurr = n;
int imin = 0;
int jmin = 0;
int node = n;
for (; node < 2 * n - 2; node++)
{
premin(d, nextp);
double dmin = 9.99e99;
for (i = 0; i >= 0; i = nextp[i])
{
if (tp[i] == fsroot)
{
continue;
}
for (j = nextp[i]; j >= 0; j = nextp[j])
{
if (tp[j] == fsroot)
{
continue;
}
double dd = dminfn(d, i, j);
if ((dd) < dmin)
{
dmin = dd;
imin = i;
jmin = j;
}
}
}

i = imin;
j = jmin;
t[tp[i]].mo = t[tp[j]].mo = node;
t[tp[i]].modist = dbranchfn(d, i, j);
t[tp[j]].modist = dbranchfn(d, j, i);
t[node].ldau = tp[i];
t[node].rdau = tp[j];
t[node].nel = t[tp[i]].nel + t[tp[j]].nel;
for (int k = 0; k >= 0; k = nextp[k])
{
tmp[k] = dnewfn(d, k, i, j, t[tp[i]].nel, t[tp[j]].nel);
}
for (int k = 0; k >= 0; k = nextp[k])
{
d[i, k] = d[k, i] = tmp[k];
}
tp[i] = node;
if (prevp[j] >= 0)
{
nextp[prevp[j]] = nextp[j];
}
if (nextp[j] >= 0)
{
prevp[nextp[j]] = prevp[j];
}
ncurr--;
}

i = 0;
j = nextp[0];
root = node;
t[tp[i]].mo = t[tp[j]].mo = t[root].mo = root;

//drootbranchfn(d, i, j, t[tp[i]].nel, t[tp[j]].nel, ref t[tp[i]].modist, ref t[tp[j]].modist);
double tpim = t[tp[i]].modist;
double tpjm = t[tp[j]].modist;
drootbranchfn(d, i, j, t[tp[i]].nel, t[tp[j]].nel, ref tpim, ref tpjm);
t[tp[i]].modist = tpim;
t[tp[j]].modist = tpjm;

t[root].ldau = tp[i];
t[root].rdau = tp[j];
t[root].modist = t[root].dep = 0.0;
t[root].nel = t[tp[i]].nel + t[tp[j]].nel;
int ntask = 0;
seqmax = depmax = 0.0;
tasklist[ntask++] = root;
while (ntask > 0)
{
i = tasklist[--ntask];
if (i >= 0)
{
t[i].dep = t[t[i].mo].dep + t[i].modist;
if (t[i].dep > depmax)
{
depmax = t[i].dep;
}
if (t[i].ldau < 0)
{
t[i].seq = seqmax++;
}
else
{
tasklist[ntask++] = -i - 1;
tasklist[ntask++] = t[i].ldau;
tasklist[ntask++] = t[i].rdau;
}
}
else
{
i = -i - 1;
t[i].seq = 0.5 * (t[t[i].ldau].seq + t[t[i].rdau].seq);
}
}
}

public int comancestor(int leafa, int leafb)
{
int i;
int j;
for (i = leafa; i != root; i = t[i].mo)
{
for (j = leafb; j != root; j = t[j].mo)
{
if (i == j)
{
break;
}
}
if (i == j)
{
break;
}
}
return i;
}

#if __UNUSED__
private void newick(Phylagglom p, MatChar str, ref string filename)
{
FILE OUT = fopen(filename, "wb");
int i;
int s;
int ntask = 0;
int n = p.n;
int root = p.root;
int[] tasklist = new int[2 * n + 1];
tasklist[ntask++] = (1 << 16) + root;
while (ntask-- > 0)
{
s = tasklist[ntask] >> 16;
i = tasklist[ntask] & 0xffff;
if (s == 1 || s == 2)
{
tasklist[ntask++] = ((s + 2) << 16) + p.t[i].mo;
if (p.t[i].ldau >= 0)
{
fprintf(OUT, "(");
tasklist[ntask++] = (2 << 16) + p.t[i].rdau;
tasklist[ntask++] = (1 << 16) + p.t[i].ldau;
}
else
{
fprintf(OUT, "%s:%f", str[i, 0], p.t[i].modist);
}
}
else if (s == 3)
{
if (ntask > 0)
{
fprintf(OUT, ",\n");
}
}
else if (s == 4)
{
if (i == root)
{
fprintf(OUT, ");\n");
}
else
{
fprintf(OUT, "):%f", p.t[i].modist);
}
}
}
fclose(OUT);
}

private void phyl2ps(ref string filename, Phylagglom ph, MatChar str, int extend, double xl, double xr, double yt, double yb)
{
int i;
int j;
double id;
double jd;
double xi;
double yi;
double xj;
double yj;
double seqmax;
double depmax;
FILE OUT = fopen(filename, "wb");
fprintf(OUT, "%%!PS\n/Courier findfont 8 scalefont setfont\n");
seqmax = ph.seqmax;
depmax = ph.depmax;
for (i = 0; i < 2 * (ph.n) - 1; i++)
{
j = ph.t[i].mo;
id = ph.t[i].dep;
jd = ph.t[j].dep;
xi = xl + (xr - xl) * id / depmax;
yi = yt - (yt - yb) * (ph.t[i].seq + 0.5) / seqmax;
xj = xl + (xr - xl) * jd / depmax;
yj = yt - (yt - yb) * (ph.t[j].seq + 0.5) / seqmax;
fprintf(OUT, "%f %f moveto %f %f lineto %f %f lineto 0 setgray stroke\n", xj, yj, xj, yi, xi, yi);
if (extend != 0)
{
if (i < ph.n)
{
fprintf(OUT, "%f %f moveto %f %f lineto 0.7 setgray stroke\n", xi, yi, xr, yi);
fprintf(OUT, "%f %f moveto (%s (%02d)) 0 setgray show\n", xr + 3.0, yi - 2.0, str[i, 0], i);
}
}
else
{
if (i < ph.n)
{
fprintf(OUT, "%f %f moveto (%s (%02d)) 0 setgray show\n", xi + 3.0, yi - 2.0, str[i, 0], i);
}
}
}
fprintf(OUT, "showpage\n\x0004");
fclose(OUT);
}
#endif
}
}

2 代码格式

using System;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
    public abstract class Phylagglom
    {
        public int n { get; set; }
        public int root { get; set; }
        public int fsroot { get; set; }
        public double seqmax { get; set; }
        public double depmax { get; set; }
        public List<Phylagglomnode> t { get; set; } = new List<Phylagglomnode>();
        public abstract void premin(double[,] d, int[] nextp);
        public abstract double dminfn(double[,] d, int i, int j);
        public abstract double dbranchfn(double[,] d, int i, int j);
        public abstract double dnewfn(double[,] d, int k, int i, int j, int ni, int nj);
        public abstract void drootbranchfn(double[,] d, int i, int j, int ni, int nj, ref double bi, ref double bj);

        public Phylagglom(double[,] dist, int fsr = -1)
        {
            this.n = dist.GetLength(0);
            this.fsroot = fsr;
            this.t = new List<Phylagglomnode>(2 * n - 1);
        }

        public void makethetree(double[,] dist)
        {
            double[,] d = Globals.CopyFrom(dist);
            int[] tp = new int[n];
            int[] nextp = new int[n];
            int[] prevp = new int[n];
            int[] tasklist = new int[2 * n + 1];
            double[] tmp = new double[n];

            int j;
            int i = 0;
            for (; i < n; i++)
            {
                nextp[i] = i + 1;
                prevp[i] = i - 1;
                tp[i] = i;
                t[i].ldau = t[i].rdau = -1;
                t[i].nel = 1;
            }

            prevp[0] = nextp[n - 1] = -1;
            int ncurr = n;
            int imin = 0;
            int jmin = 0;
            int node = n;
            for (; node < 2 * n - 2; node++)
            {
                premin(d, nextp);
                double dmin = 9.99e99;
                for (i = 0; i >= 0; i = nextp[i])
                {
                    if (tp[i] == fsroot)
                    {
                        continue;
                    }
                    for (j = nextp[i]; j >= 0; j = nextp[j])
                    {
                        if (tp[j] == fsroot)
                        {
                            continue;
                        }
                        double dd = dminfn(d, i, j);
                        if ((dd) < dmin)
                        {
                            dmin = dd;
                            imin = i;
                            jmin = j;
                        }
                    }
                }

                i = imin;
                j = jmin;
                t[tp[i]].mo = t[tp[j]].mo = node;
                t[tp[i]].modist = dbranchfn(d, i, j);
                t[tp[j]].modist = dbranchfn(d, j, i);
                t[node].ldau = tp[i];
                t[node].rdau = tp[j];
                t[node].nel = t[tp[i]].nel + t[tp[j]].nel;
                for (int k = 0; k >= 0; k = nextp[k])
                {
                    tmp[k] = dnewfn(d, k, i, j, t[tp[i]].nel, t[tp[j]].nel);
                }
                for (int k = 0; k >= 0; k = nextp[k])
                {
                    d[i, k] = d[k, i] = tmp[k];
                }
                tp[i] = node;
                if (prevp[j] >= 0)
                {
                    nextp[prevp[j]] = nextp[j];
                }
                if (nextp[j] >= 0)
                {
                    prevp[nextp[j]] = prevp[j];
                }
                ncurr--;
            }

            i = 0;
            j = nextp[0];
            root = node;
            t[tp[i]].mo = t[tp[j]].mo = t[root].mo = root;

            //drootbranchfn(d, i, j, t[tp[i]].nel, t[tp[j]].nel, ref t[tp[i]].modist, ref t[tp[j]].modist);
            double tpim = t[tp[i]].modist;
            double tpjm = t[tp[j]].modist;
            drootbranchfn(d, i, j, t[tp[i]].nel, t[tp[j]].nel, ref tpim, ref tpjm);
            t[tp[i]].modist = tpim;
            t[tp[j]].modist = tpjm;

            t[root].ldau = tp[i];
            t[root].rdau = tp[j];
            t[root].modist = t[root].dep = 0.0;
            t[root].nel = t[tp[i]].nel + t[tp[j]].nel;
            int ntask = 0;
            seqmax = depmax = 0.0;
            tasklist[ntask++] = root;
            while (ntask > 0)
            {
                i = tasklist[--ntask];
                if (i >= 0)
                {
                    t[i].dep = t[t[i].mo].dep + t[i].modist;
                    if (t[i].dep > depmax)
                    {
                        depmax = t[i].dep;
                    }
                    if (t[i].ldau < 0)
                    {
                        t[i].seq = seqmax++;
                    }
                    else
                    {
                        tasklist[ntask++] = -i - 1;
                        tasklist[ntask++] = t[i].ldau;
                        tasklist[ntask++] = t[i].rdau;
                    }
                }
                else
                {
                    i = -i - 1;
                    t[i].seq = 0.5 * (t[t[i].ldau].seq + t[t[i].rdau].seq);
                }
            }
        }

        public int comancestor(int leafa, int leafb)
        {
            int i;
            int j;
            for (i = leafa; i != root; i = t[i].mo)
            {
                for (j = leafb; j != root; j = t[j].mo)
                {
                    if (i == j)
                    {
                        break;
                    }
                }
                if (i == j)
                {
                    break;
                }
            }
            return i;
        }

#if __UNUSED__
        private void newick(Phylagglom p, MatChar str, ref string filename)
        {
            FILE OUT = fopen(filename, "wb");
            int i;
            int s;
            int ntask = 0;
            int n = p.n;
            int root = p.root;
            int[] tasklist = new int[2 * n + 1];
            tasklist[ntask++] = (1 << 16) + root;
            while (ntask-- > 0)
            {
                s = tasklist[ntask] >> 16;
                i = tasklist[ntask] & 0xffff;
                if (s == 1 || s == 2)
                {
                    tasklist[ntask++] = ((s + 2) << 16) + p.t[i].mo;
                    if (p.t[i].ldau >= 0)
                    {
                        fprintf(OUT, "(");
                        tasklist[ntask++] = (2 << 16) + p.t[i].rdau;
                        tasklist[ntask++] = (1 << 16) + p.t[i].ldau;
                    }
                    else
                    {
                        fprintf(OUT, "%s:%f", str[i, 0], p.t[i].modist);
                    }
                }
                else if (s == 3)
                {
                    if (ntask > 0)
                    {
                        fprintf(OUT, ",\n");
                    }
                }
                else if (s == 4)
                {
                    if (i == root)
                    {
                        fprintf(OUT, ");\n");
                    }
                    else
                    {
                        fprintf(OUT, "):%f", p.t[i].modist);
                    }
                }
            }
            fclose(OUT);
        }

        private void phyl2ps(ref string filename, Phylagglom ph, MatChar str, int extend, double xl, double xr, double yt, double yb)
        {
            int i;
            int j;
            double id;
            double jd;
            double xi;
            double yi;
            double xj;
            double yj;
            double seqmax;
            double depmax;
            FILE OUT = fopen(filename, "wb");
            fprintf(OUT, "%%!PS\n/Courier findfont 8 scalefont setfont\n");
            seqmax = ph.seqmax;
            depmax = ph.depmax;
            for (i = 0; i < 2 * (ph.n) - 1; i++)
            {
                j = ph.t[i].mo;
                id = ph.t[i].dep;
                jd = ph.t[j].dep;
                xi = xl + (xr - xl) * id / depmax;
                yi = yt - (yt - yb) * (ph.t[i].seq + 0.5) / seqmax;
                xj = xl + (xr - xl) * jd / depmax;
                yj = yt - (yt - yb) * (ph.t[j].seq + 0.5) / seqmax;
                fprintf(OUT, "%f %f moveto %f %f lineto %f %f lineto 0 setgray stroke\n", xj, yj, xj, yi, xi, yi);
                if (extend != 0)
                {
                    if (i < ph.n)
                    {
                        fprintf(OUT, "%f %f moveto %f %f lineto 0.7 setgray stroke\n", xi, yi, xr, yi);
                        fprintf(OUT, "%f %f moveto (%s (%02d)) 0 setgray show\n", xr + 3.0, yi - 2.0, str[i, 0], i);
                    }
                }
                else
                {
                    if (i < ph.n)
                    {
                        fprintf(OUT, "%f %f moveto (%s (%02d)) 0 setgray show\n", xi + 3.0, yi - 2.0, str[i, 0], i);
                    }
                }
            }
            fprintf(OUT, "showpage\n\x0004");
            fclose(OUT);
        }
#endif
    }
}