SHAPPE/sampleTargetDistance.m at main · Craig-Lab/SHAPPE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
function s = sampleTargetDistance(Neigh,distributionsObs,distributionsRandom,copulaInformation,types,...
    samplingindexes,types_names,correlatedSamplingIdx,uncorrelatedSamplingIdx,exclusionDistance)

    copulaType = lower(string(copulaInformation.kernel));
    rho = copulaInformation.rho;
    if(copulaType == "t")
        nu = copulaInformation.nu;
    end

    Ntypes = length(types);
    s = zeros(1,Ntypes*Neigh);

    % sample a copula row
    if(copulaType=="t")
        U = copularnd('t',rho,nu,1);
    elseif(copulaType=="gaussian")
        U = copularnd('Gaussian',rho,1);
    else
        % No copula: sample indpendant probabilities
        U = unifrnd(0,1,[1 size(rho,2)]);
    end

    for nD=correlatedSamplingIdx
        % find which type this nD belongs to
        current_dist = [];
        for t = 1:Ntypes
            ct_name  = types_names(t);

            if nD >= samplingindexes(1,t) && nD <= samplingindexes(2,t)
                if(isfield(distributionsObs,ct_name))
                    current_dist = distributionsObs.(ct_name).distribution{nD - samplingindexes(1,t) + 1};
                    cumProb = distributionsObs.(ct_name).cum_probs{nD - samplingindexes(1,t) + 1};

                    if isempty(current_dist) || ~any(current_dist > 0)
                        current_dist = [];
                        break
                    end

                    if nD == samplingindexes(1,t)
                        cumulDistance = 0;
                    end
                end
                break
            end
        end

        if isempty(current_dist)
            continue
        end

        ind_U = find(correlatedSamplingIdx==nD, 1, "first");
        copulaInformation.sampleIndex = ind_U;

        foundDistance = false;
        while ~foundDistance
            if numel(correlatedSamplingIdx) > 1
                u = U(1, ind_U);
            else
                u = rand();
            end
            u = min(max(u, realmin), 1-realmin);
            d = getDistanceProbability(current_dist, cumProb,u);
            foundDistance = isfinite(d) && (d >= 0);
        end

        if nD < (Ntypes*Neigh + 1)
            cumulDistance = cumulDistance + d;
            d = cumulDistance;
        end

        s(1,nD) = d;
    end % End of for nD=correlatedSamplingIdx

    for nD=uncorrelatedSamplingIdx
        current_dist = [];

        for t = 1:Ntypes
            ct_name  = types_names(t);

            if nD >= samplingindexes(1,t) && nD <= samplingindexes(2,t)
                current_dist = distributionsRandom.(ct_name).distribution{nD - samplingindexes(1,t) + 1};
                cumProb = distributionsRandom.(ct_name).cum_probs{nD - samplingindexes(1,t) + 1};

                if isempty(current_dist) || ~any(current_dist > 0)
                    current_dist = [];
                    break
                end

                if nD == samplingindexes(1,t)
                    cumulDistance = 0;
                end
                break
            end
        end

        if isempty(current_dist)
            continue
        end

        foundDistance = false;
        while ~foundDistance
            u = rand();
            u = min(max(u, realmin), 1-realmin);
            d = getDistanceProbability(current_dist, cumProb,u);
            foundDistance = isfinite(d) && (d >= 0);
        end

        cumulDistance = cumulDistance + d;
        s(1,nD) = cumulDistance;
    end % End of for nD=uncorrelatedSamplingIdx

end