Ignore:
Timestamp:
12/08/20 08:45:53 (4 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 25834

Location:
issm/trunk
Files:
11 edited
3 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/m/classes/qmu/continuous_design.m

    r24686 r25836  
    220220                        pairs_per_variable=[];
    221221        end % }}}
     222                function checkconsistency(self,md,solution,analyses) % {{{
     223                end % }}}
    222224
    223225        end

  • issm/trunk/src/m/classes/qmu/histogram_bin_uncertain.m

    r24686 r25836  
    1 % %  definition for the histogram_bin_uncertain class.
     1%HISTOGRAM BIN UNCERTAIN class definition
    22%
    3 %  [hbu]=histogram_bin_uncertain(varargin)
     3%       Usage:
     4%   hbu=histogram_bin_uncertain('descriptor',descriptor,'pairs_per_variable',pairs_per_variable,...
     5%           'abscissas',abscissas,'counts',counts,'partition',partition)
     6%     
     7%      where hbu is the histogram_bin_uncertain object returned by the constructor, pairs_per_variable the
     8%      size of the histogram, 'abscissas' and 'counts' describe the histogram itself.
     9%      If the variable is distributed, then a partition vector can be supplied, which can be a
     10%      partition vector over elements or vertices. In which case counts,
    411%
    5 %  where the required varargin are:
    6 %    descriptor    (char, description, '')
    7 %    pairs_per_variable          (double vector, [])
    8 %    abscissas          (double vector, [])
    9 %    counts          (int vector, [])
    10 %
    11 %  note that zero arguments constructs a default instance; one
    12 %  argument of the class copies the instance; and three or more
    13 %  arguments constructs a new instance from the arguments.
    14 %
     12%   Example:
     13%      md.qmu.variables.giaid=histogram_bin_uncertain('descriptor','GiammeModelId','pairs_per_variable',3,...
     14%           'count',[.6 .4 0],''abscissas',[1 2 3]);
     15%      md.qmu.variables.surfaceloadid=histogram_bin_uncertain(...
     16%           'descriptor','distributed_SurfaceloadModelId','pairs_per_variable',[2 3 4],...
     17%           'counts',{[1 0],[.6 .4 0],[.4 .4 .2 0]},...
     18%           'abscissas',{[1 2],[1 2 3],[1 2 3 4]},...
     19%           'partition',partition_vector);
     20
    1521classdef histogram_bin_uncertain
    16     properties
    17         descriptor='';
    18                 pairs_per_variable=[];
    19         abscissas = [];
    20         counts = [];
    21     end
     22        properties
     23                descriptor='';
     24                pairs_per_variable={};
     25                abscissas = {};
     26                counts = [];
     27                partition=[];
     28        end
     29        methods
     30                function self=histogram_bin_uncertain(varargin) % {{{
     31               
     32                        %recover options:
     33                        options=pairoptions(varargin{:});
    2234
    23     methods
    24         function [hbu]=histogram_bin_uncertain(varargin) % {{{
     35                        %initialize fields:
     36                        self.descriptor=getfieldvalue(options,'descriptor');
     37                        self.pairs_per_variable=getfieldvalue(options,'pairs_per_variable');
     38                        self.abscissas=getfieldvalue(options,'abscissas');
     39                        self.counts=getfieldvalue(options,'counts');
    2540
    26             switch nargin
    27                 case 0 %  create a default object
    28                 case 1 %  copy the object
    29                     if isa(varargin{1},'histogram_bin_uncertain')
    30                         hbu=varargin{1};
    31                     else
    32                         error('Object ''%s'' is a ''%s'' class object, not ''%s''.',...
    33                             inputname(1),class(varargin{1}),'histogram_bin_uncertain');
    34                     end
    35                 case {2,3} %  not enough arguments
    36                     error('Construction of ''%s'' class object requires at least %d inputs.',...
    37                         'histogram_bin_uncertain',4)
    38                 case 4 %
    39                                         %  create the object from the input
    40                                         hbu = histogram_bin_uncertain;
    41                                         hbu.descriptor=varargin{1};
    42                                         hbu.pairs_per_variable=varargin{2};
    43                                         hbu.abscissas=varargin{3};
    44                                         hbu.counts=varargin{4};
    45 
    46                 otherwise
    47                                         error('Construction of histogram_bin_uncertain class object requires three arguments, descriptor, abscissas and counts');
    48             end
    49 
    50         end % }}}
    51         function []=disp(hbu) % {{{
    52 
    53 %  display the object
    54 
    55             disp(sprintf('\n'));
    56             for i=1:numel(hbu)
    57                 disp(sprintf('class ''%s'' object ''%s%s'' = \n',...
    58                     class(hbu),inputname(1),string_dim(hbu,i)));
    59                 disp(sprintf('    descriptor: ''%s'''  ,hbu(i).descriptor));
    60                 disp(sprintf('          pairs_per_variable: %g'      ,hbu(i).pairs_per_variable));
    61                 disp(sprintf('          abscissas: %g'      ,hbu(i).abscissas));
    62                 disp(sprintf('        counts: %g'      ,hbu(i).counts));
    63             end
    64 
    65         end % }}}
    66         function [desc]  =prop_desc(hbu,dstr) % {{{
     41                        %if the variable is distributed,  a partition vector should have been
     42                        %supplied, and that partition vector should have as many partitions
     43                        %as the pairs_per_variable, abscissas and counts arrays:
     44                        if self.isdistributed(),
     45                                self.partition=getfieldvalue(options,'partition');
     46                                npart=qmupart2npart(self.partition);
     47                                if npart~=length(self.pairs_per_variable),
     48                                        error(sprintf('histogram_bin_uncertain constructor: for the distributed variable  %s the number of pairs_per_variable arrays (%i) should be equal to the number of partitions (%i)',self.descriptor,length(self.pairs_per_variable),npart));
     49                                end
     50                                if npart~=length(self.abscissas),
     51                                        error(sprintf('histogram_bin_uncertain constructor: for the distributed variable  %s the number of abscissas arrays (%i) should be equal to the number of partitions (%i)',self.descriptor,length(self.abscissas),npart));
     52                                end
     53                                if npart~=length(self.counts),
     54                                        error(sprintf('histogram_bin_uncertain constructor: for the distributed variable  %s the size of counts (%i) should be equal to the number of partitions (%i)',self.descriptor,length(self.counts),npart));
     55                                end
     56                        end
     57                end % }}}
     58                function md=checkconsistency(self,md,solution,analyses) % {{{
     59                end % }}}
     60                        function disp(self) % {{{
     61                        disp(sprintf('   histogram_bin uncertain variable: '));
     62                        fielddisplay(self,'descriptor','name tag');
     63                        fielddisplay(self,'pairs_per_variable','number of bins in histogram');
     64                        fielddisplay(self,'abscissas','abscissas for histogram');
     65                        fielddisplay(self,'counts','probabilities for histogram');
     66                        if ~isempty(self.partition),
     67                                fielddisplay(self,'partition','partition vector defining where sampling will occur');
     68                        end
     69                end
     70                %}}}
     71                %virtual functions needed by qmu processing algorithms:
     72                %implemented:
     73        function [desc]=prop_desc(hbu,dstr) % {{{
    6774            desc=cell(1,numel(hbu));
    6875            for i=1:numel(hbu)
    … …  
    7986            desc=allempty(desc);
    8087        end  %}}}
    81         function [initpt]=prop_initpt(hbu) % {{{
    82             initpt=[];
    83         end % }}}
    84         function [lower] =prop_lower(hbu) % {{{
    85             lower=[];
    86         end % }}}
    87         function [upper] =prop_upper(hbu) % {{{
    88             upper=[];
    89         end % }}}
    90         function [mean]  =prop_mean(hbu) % {{{
     88        function [mean]=prop_mean(hbu) % {{{
    9189            mean=[];
    9290        end % }}}
    … …  
    9492            stddev=[];
    9593        end % }}}
     94        function [lower]=prop_lower(hbu) % {{{
     95            lower=[];
     96        end % }}}
     97        function [upper]=prop_upper(hbu) % {{{
     98            upper=[];
     99        end % }}}
     100                function [pairs_per_variable] =prop_pairs_per_variable(hbu) % {{{
     101                        pairs_per_variable=[];
     102            for i=1:numel(hbu)
     103                pairs_per_variable=[pairs_per_variable hbu(i).pairs_per_variable];
     104            end
     105        end % }}}
     106                function [abscissas]=prop_abscissas(hbu) % {{{
     107                        abscissas=[];
     108                        for i=1:numel(hbu)
     109                                abscissas=[abscissas hbu(i).abscissas'];
     110                        end
     111                        abscissas=allequal(abscissas,-Inf);
     112                end % }}}
     113                function [counts] =prop_counts(hbu) % {{{
     114                        counts=[];
     115                        for i=1:numel(hbu)
     116                                counts=[counts hbu(i).counts'];
     117                        end
     118                        counts=allequal(counts,-Inf);
     119        end % }}}
     120        function [initpt]=prop_initpt(hbu) % {{{
     121            initpt=[];
     122        end % }}}
    96123        function [initst]=prop_initst(hbu) % {{{
    97124            initst=[];
    98125        end % }}}
    99         function [stype] =prop_stype(hbu) % {{{
     126        function [stype]=prop_stype(hbu) % {{{
    100127            stype={};
    101128        end % }}}
    102         function [scale] =prop_scale(hbu) % {{{
     129        function [scale]=prop_scale(hbu) % {{{
    103130            scale=[];
    104131        end % }}}
    105                 function [abscissas] =prop_abscissas(hbu) % {{{
    106                 abscissas=[];
    107                 for i=1:numel(hbu)
    108                         abscissas=[abscissas hbu(i).abscissas];
    109                 end
    110                 abscissas=allequal(abscissas,-Inf);
    111 
    112         end % }}}
    113                 function [pairs_per_variable] =prop_pairs_per_variable(hbu) % {{{
    114                         pairs_per_variable=zeros(1,numel(hbu));
    115             for i=1:numel(hbu)
    116                 pairs_per_variable(i)=hbu(i).pairs_per_variable;
    117             end
    118             pairs_per_variable=allequal(pairs_per_variable,-Inf);
    119         end % }}}
    120                 function [counts] =prop_counts(hbu) % {{{
    121                 counts=[];
    122                 for i=1:numel(hbu)
    123                         counts=[counts hbu(i).counts];
    124                 end
    125                 counts=allequal(counts,-Inf);
    126 
    127         end % }}}
     132                function scaled=isscaled(self) % {{{
     133                        if strncmp(self.descriptor,'scaled_',7),
     134                                scaled=1;
     135                        else
     136                                scaled=0;
     137                        end
     138                end % }}}
     139                function distributed=isdistributed(self) % {{{
     140                        if strncmp(self.descriptor,'distributed_',12),
     141                                distributed=1;
     142                        else
     143                                distributed=0;
     144                        end
     145                end % }}}
    128146        end
    129147    methods (Static)
    130148        function []=dakota_write(fidi,dvar) % {{{
     149                        % collect only the variables of the appropriate class
     150                        hbu=struc_class(dvar,'histogram_bin_uncertain');
    131151
    132 %  collect only the variables of the appropriate class
    133 
    134             hbu=struc_class(dvar,'histogram_bin_uncertain');
    135 
    136 %  write variables
    137 
     152                        % write variables
    138153            vlist_write(fidi,'histogram_bin_uncertain','hbu',hbu);
    139154        end % }}}

  • issm/trunk/src/m/classes/qmu/linear_inequality_constraint.py

    r24313 r25836  
    1313  where the required args are:
    1414    matrix        (double row, variable coefficients, float('NaN'))
    15     lower         (double vector, lower bounds, -np.Inf)
     15    lower         (double vector, lower bounds, -np.inf)
    1616    upper         (double vector, upper bounds, 0.)
    1717  and the optional args and defaults are:
    … …  
    2525    def __init__(self):
    2626        self.matrix = np.array([[float('NaN')]])
    27         self.lower = -np.Inf
     27        self.lower = -np.inf
    2828        self.upper = 0.
    2929        self.scale_type = 'none'
    … …  
    131131            lower[i] = lic[i].lower
    132132
    133         lower = allequal(lower, -np.Inf)
     133        lower = allequal(lower, -np.inf)
    134134
    135135        return lower

  • issm/trunk/src/m/classes/qmu/normal_uncertain.m

    r24686 r25836  
     1%NORMAL_UNCERTAIN class definition
    12%
    2 %  definition for the normal_uncertain class.
     3%   Usage:
     4%      nuv=normal_uncertain('descriptor',descriptor,'mean',mean,'stddev',stddev,'partition',partition);
     5%      where nuv is the normal_uncertain object returned by the constructor, mean and stddev are self
     6%      explanatory.  partition is the partition vector for distributed variables. Can be a partition
     7%      vector over elements or vertices.
    38%
    4 %  [nuv]=normal_uncertain(varargin)
     9%   Example:
     10%      md.qmu.variables.rheology=normal_uncertain('descriptor','RheologyBBar','mean',1,'stddev',.05);
     11%      md.qmu.variables.rheology=normal_uncertain('descriptor','scaled_RheologyBBar','mean',1,'stddev',.05,'partition',vpartition);
    512%
    6 %  where the required varargin are:
    7 %    descriptor    (char, description, '')
    8 %    mean          (double, mean, NaN)
    9 %    stddev        (double, standard deviation, NaN)
    10 %  and the optional varargin and defaults are:
    11 %    lower         (double, lower bound, -Inf)
    12 %    upper         (double, upper bound,  Inf)
    13 %
    14 %  note that zero arguments constructs a default instance; one
    15 %  argument of the class copies the instance; and three or more
    16 %  arguments constructs a new instance from the arguments.
    17 %
    18 %  "Copyright 2009, by the California Institute of Technology.
    19 %  ALL RIGHTS RESERVED. United States Government Sponsorship
    20 %  acknowledged. Any commercial use must be negotiated with
    21 %  the Office of Technology Transfer at the California Institute
    22 %  of Technology.  (J. Schiermeier, NTR 47078)
    23 %
    24 %  This software may be subject to U.S. export control laws.
    25 %  By accepting this  software, the user agrees to comply with
    26 %  all applicable U.S. export laws and regulations. User has the
    27 %  responsibility to obtain export licenses, or other export
    28 %  authority as may be required before exporting such information
    29 %  to foreign countries or providing access to foreign persons."
    30 %
     13
    3114classdef normal_uncertain
    32     properties
    33         descriptor='';
    34         mean      = NaN;
    35         stddev    = NaN;
    36         lower     =-Inf;
    37         upper     = Inf;
    38     end
     15        properties
     16                descriptor      = '';
     17                mean            = NaN;
     18                stddev          = NaN;
     19                partition       = [];
     20                nsteps          = 0;
     21        end
     22        methods
     23                function self=normal_uncertain(varargin) %constructor {{{
    3924
    40     methods
    41         function [nuv]=normal_uncertain(varargin)
     25                        %recover options:
     26                        options=pairoptions(varargin{:});
    4227
    43             switch nargin
     28                        %initialize fields:
     29                        self.descriptor=getfieldvalue(options,'descriptor');
     30                        self.mean=getfieldvalue(options,'mean');
     31                        self.stddev=getfieldvalue(options,'stddev');
    4432
    45 %  create a default object
     33                        %if the variable is scaled,  a partition vector should have been
     34                        %supplied, and that partition vector should have as many partitions
     35                        %as the mean and stddev vectors:
     36                        if self.isscaled(),
     37                                self.partition=getfieldvalue(options,'partition');
     38                                self.nsteps=getfieldvalue(options,'nsteps',1);
     39                                npart=qmupart2npart(self.partition);
     40                                if npart~=size(self.mean,1),
     41                                        error(['normal_uncertain constructor: for the scaled variable ' self.descriptor ' the row size of the mean field should be identical to the number of partitions']);
     42                                end
     43                                if npart~=size(self.stddev,1),
     44                                        error(['normal_uncertain constructor: for the scaled variable ' self.descriptor ' the row size of the stddev field should be identical to the number of partitions']);
     45                                end
     46                                if self.nsteps~=size(self.mean,2),
     47                                        error(['normal_uncertain constructor: for the scaled variable ' self.descriptor ' the col size of the mean field should be identical to the number of time steps']);
     48                                end
     49                                if self.nsteps~=size(self.stddev,2),
     50                                        error(['normal_uncertain constructor: for the scaled variable ' self.descriptor ' the col size of the stddev field should be identical to the number of time steps']);
     51                                end
    4652
    47                 case 0
     53                        end
    4854
    49 %  copy the object
     55                end %}}}
     56                function disp(self) % {{{
     57                        disp(sprintf('   normal uncertain variable: '));
     58                        fielddisplay(self,'descriptor','name tag');
     59                        fielddisplay(self,'mean','pdf mean');
     60                        fielddisplay(self,'stddev','pdf standard deviation');
     61                        if ~isempty(self.partition),
     62                                fielddisplay(self,'partition','partition vector defining where sampling will occur');
     63                        end
     64                        fielddisplay(self,'nsteps','number of time steps');
     65                end
     66                %}}}
     67                function md=checkconsistency(self,md,solution,analyses) % {{{
    5068
    51                 case 1
    52                     if isa(varargin{1},'normal_uncertain')
    53                         nuv=varargin{1};
    54                     else
    55                         error('Object ''%s'' is a ''%s'' class object, not ''%s''.',...
    56                             inputname(1),class(varargin{1}),'normal_uncertain');
    57                     end
     69                        md = checkfield(md,'field',self.mean,'fieldname','normal_uncertain.mean','NaN',1,'Inf',1,'>=',0);
     70                        md = checkfield(md,'field',self.stddev,'fieldname','normal_uncertain.stddev','NaN',1,'Inf',1,'>=',0);
     71                        if self.isscaled(),
     72                                if isempty(self.partition),
     73                                        error('normal_uncertain is a scaled variable, but it''s missing a partition vector');
     74                                end
     75                                %better have a partition vector that has as many partitions as stddev's size:
     76                                if size(self.stddev,1)~=partition_npart(self.partition),
     77                                        error('normal_uncertain error message: row size of stddev and partition size should be identical');
     78                                end
     79                                if size(self.mean,1)~=partition_npart(self.partition),
     80                                        error('normal_uncertain error message: row size of mean and partition size should be identical');
     81                                end
     82                                %we need as steps in stddev and mean as there are time steps:
     83                                if size(self.stddev,2)~=self.nsteps,
     84                                        error('normal_uncertain error message: col size of stddev and number of time steps should be identical');
     85                                end
     86                                if size(self.mean,2)~=self.nsteps,
     87                                        error('normal_uncertain error message: col size of mean and number of time steps should be identical');
     88                                end
    5889
    59 %  not enough arguments
    60 
    61                 case 2
    62                     error('Construction of ''%s'' class object requires at least %d inputs.',...
    63                         'normal_uncertain',3)
    64 
    65 %  create the object from the input
    66 
    67                 otherwise
    68                     asizec=num2cell(array_size(varargin{1:min(nargin,5)}));
    69                     nuv(asizec{:})=normal_uncertain;
    70                     clear asizec
    71 
    72                     if ischar(varargin{1})
    73                         varargin{1}=cellstr(varargin{1});
    74                     end
    75                     for i=1:numel(nuv)
    76                         if (numel(varargin{1}) > 1)
    77                             nuv(i).descriptor=varargin{1}{i};
    78                         else
    79                             if numel(nuv)==1,
    80                                                                 nuv(i).descriptor=char(varargin{1});
    81                                                         else
    82                                                                 nuv(i).descriptor=[char(varargin{1}) num2str(i)];
    83                                                         end
    84                         end
    85                         if (numel(varargin{2}) > 1)
    86                             nuv(i).mean      =varargin{2}(i);
    87                         else
    88                             nuv(i).mean      =varargin{2};
    89                         end
    90                         if (numel(varargin{3}) > 1)
    91                             nuv(i).stddev    =varargin{3}(i);
    92                         else
    93                             nuv(i).stddev    =varargin{3};
    94                         end
    95                     end
    96 
    97                     if (nargin >= 4)
    98                         for i=1:numel(nuv)
    99                             if (numel(varargin{4}) > 1)
    100                                 nuv(i).lower     =varargin{4}(i);
    101                             else
    102                                 nuv(i).lower     =varargin{4};
    103                             end
    104                         end
    105                         if (nargin >= 5)
    106                             for i=1:numel(nuv)
    107                                 if (numel(varargin{5}) > 1)
    108                                     nuv(i).upper     =varargin{5}(i);
    109                                 else
    110                                     nuv(i).upper     =varargin{5};
    111                                 end
    112                             end
    113                             if (nargin > 5)
    114                                 warning('normal_uncertain:extra_arg',...
    115                                     'Extra arguments for object of class ''%s''.',...
    116                                     class(nuv));
    117                             end
    118                         end
    119                     end
    120             end
    121 
    122         end
    123 
    124         function []=disp(nuv)
    125 
    126 %  display the object
    127 
    128             disp(sprintf('\n'));
    129             for i=1:numel(nuv)
    130                 disp(sprintf('class ''%s'' object ''%s%s'' = \n',...
    131                     class(nuv),inputname(1),string_dim(nuv,i)));
    132                 disp(sprintf('    descriptor: ''%s'''  ,nuv(i).descriptor));
    133                 disp(sprintf('          mean: %g'      ,nuv(i).mean));
    134                 disp(sprintf('        stddev: %g'      ,nuv(i).stddev));
    135                 disp(sprintf('         lower: %g'      ,nuv(i).lower));
    136                 disp(sprintf('         upper: %g\n'    ,nuv(i).upper));
    137             end
    138 
    139         end
    140 
    141         function [desc]  =prop_desc(nuv,dstr)
    142             desc=cell(1,numel(nuv));
    143             for i=1:numel(nuv)
    144                 if ~isempty(nuv(i).descriptor)
    145                     desc(i)=cellstr(nuv(i).descriptor);
    146                 elseif ~isempty(inputname(1))
    147                     desc(i)=cellstr([inputname(1) string_dim(nuv,i,'vector')]);
    148                 elseif exist('dstr','var')
    149                     desc(i)=cellstr([dstr         string_dim(nuv,i,'vector')]);
    150                 else
    151                     desc(i)=cellstr(['nuv'        string_dim(nuv,i,'vector')]);
    152                 end
    153             end
    154             desc=allempty(desc);
    155         end
    156         function [initpt]=prop_initpt(nuv)
    157             initpt=[];
    158         end
    159         function [lower] =prop_lower(nuv)
    160             lower=zeros(1,numel(nuv));
    161             for i=1:numel(nuv)
    162                 lower(i)=nuv(i).lower;
    163             end
    164             lower=allequal(lower,-Inf);
    165         end
    166         function [upper] =prop_upper(nuv)
    167             upper=zeros(1,numel(nuv));
    168             for i=1:numel(nuv)
    169                 upper(i)=nuv(i).upper;
    170             end
    171             upper=allequal(upper, Inf);
    172         end
    173         function [mean]  =prop_mean(nuv)
    174             mean=zeros(1,numel(nuv));
    175             for i=1:numel(nuv)
    176                 mean(i)=nuv(i).mean;
    177             end
    178         end
    179         function [stddev]=prop_stddev(nuv)
    180             stddev=zeros(1,numel(nuv));
    181             for i=1:numel(nuv)
    182                 stddev(i)=nuv(i).stddev;
    183             end
    184         end
    185         function [initst]=prop_initst(nuv)
    186             initst=[];
    187         end
    188         function [stype] =prop_stype(nuv)
    189             stype={};
    190         end
    191         function [scale] =prop_scale(nuv)
    192             scale=[];
    193         end
    194                 function [abscissas] =prop_abscissas(hbu) % {{{
    195             abscissas=[];
    196         end % }}}
    197         function [counts] =prop_counts(hbu) % {{{
    198             counts=[];
    199         end % }}}
    200         function [pairs_per_variable] =prop_pairs_per_variable(hbu) % {{{
     90                                md = checkfield(md,'field',self.partition,'fieldname','normal_uncertain.partition','NaN',1,'Inf',1,'>=',-1,'numel',[md.mesh.numberofvertices,md.mesh.numberofelements]);
     91                                if size(self.partition,2)>1,
     92                                        error('normal_uncertain error message: partition should be a column vector');
     93                                end
     94                                partcheck=unique(self.partition);
     95                                partmin=min(partcheck);
     96                                partmax=max(partcheck);
     97                                if partmax<-1,
     98                                        error('normal_uncertain error message: partition vector''s min value should be -1 (for no partition), or start at 0');
     99                                end
     100                                nmax=max(md.mesh.numberofelements,md.mesh.numberofvertices);
     101                                if partmax>nmax,
     102                                        error('normal_uncertain error message: partition vector''s values cannot go over the number of vertices or elements');
     103                                end
     104                        end
     105                end % }}}
     106                %virtual functions needed by qmu processing algorithms:
     107                %implemented:
     108                function [desc]=prop_desc(nuv,dstr) % {{{
     109                        desc=cell(1,numel(nuv));
     110                        for i=1:numel(nuv)
     111                                if ~isempty(nuv(i).descriptor)
     112                                        desc(i)=cellstr(nuv(i).descriptor);
     113                                elseif ~isempty(inputname(1))
     114                                        desc(i)=cellstr([inputname(1) string_dim(nuv,i,'vector')]);
     115                                elseif exist('dstr','var')
     116                                        desc(i)=cellstr([dstr         string_dim(nuv,i,'vector')]);
     117                                else
     118                                        desc(i)=cellstr(['nuv'        string_dim(nuv,i,'vector')]);
     119                                end
     120                        end
     121                        desc=allempty(desc);
     122                end %}}}
     123                function [mean]=prop_mean(nuv) % {{{
     124                        mean=zeros(1,numel(nuv));
     125                        for i=1:numel(nuv)
     126                                mean(i)=nuv(i).mean;
     127                        end
     128                end % }}}
     129                function [stddev]=prop_stddev(nuv) % {{{
     130                        stddev=zeros(1,numel(nuv));
     131                        for i=1:numel(nuv)
     132                                stddev(i)=nuv(i).stddev;
     133                        end
     134                end % }}}
     135                function [lower]=prop_lower(nuv) % {{{
     136                        lower=[];
     137                end % }}}
     138                function [upper]=prop_upper(nuv) % {{{
     139                        upper=[];
     140                end % }}}
     141                function [abscissas]=prop_abscissas(hbu) % {{{
     142                        abscissas=[];
     143                end % }}}
     144                function [counts]=prop_counts(hbu) % {{{
     145                        counts=[];
     146                end % }}}
     147                function [pairs_per_variable]=prop_pairs_per_variable(hbu) % {{{
    201148                        pairs_per_variable=[];
    202         end % }}}
    203 
    204     end
    205 
    206     methods (Static)
    207         function []=dakota_write(fidi,dvar)
    208 
    209 %  collect only the variables of the appropriate class
    210 
    211             nuv=struc_class(dvar,'normal_uncertain');
    212 
    213 %  write variables
    214 
    215             vlist_write(fidi,'normal_uncertain','nuv',nuv);
    216         end
    217     end
     149                end % }}}
     150                function [initpt]=prop_initpt(nuv) % {{{
     151                        initpt=[];
     152                end % }}}
     153                function [initst]=prop_initst(nuv) % {{{
     154                        initst=[];
     155                end % }}}
     156                function [stype]=prop_stype(nuv) % {{{
     157                        stype={};
     158                end % }}}
     159                function [scale]=prop_scale(nuv) % {{{
     160                        scale=[];
     161                end % }}}
     162                %new methods:
     163                function distributed=isdistributed(self) % {{{
     164                        if strncmp(self.descriptor,'distributed_',12),
     165                                distributed=1;
     166                        else
     167                                distributed=0;
     168                        end
     169                end % }}}
     170                function scaled=isscaled(self) % {{{
     171                        if strncmp(self.descriptor,'scaled_',7),
     172                                scaled=1;
     173                        else
     174                                scaled=0;
     175                        end
     176                end % }}}
     177        end
     178        methods (Static)
     179                function []=dakota_write(fidi,dvar) % {{{
     180                        %  collect only the variables of the appropriate class
     181                        nuv=struc_class(dvar,'normal_uncertain');
     182                        %  write variables
     183                        vlist_write(fidi,'normal_uncertain','nuv',nuv);
     184                end % }}}
     185        end
    218186end

  • issm/trunk/src/m/classes/qmu/normal_uncertain.py

    r24313 r25836  
    11import numpy as np
     2
    23from MatlabArray import *
     4from MatlabFuncs import *
     5from fielddisplay import fielddisplay
     6from pairoptions import pairoptions
     7from partition_npart import *
     8from qmupart2npart import qmupart2npart
    39
    410
    511class normal_uncertain(object):
    6     '''
    7   definition for the normal_uncertain class.
    8 
    9   [nuv] = normal_uncertain.normal_uncertain(args)
    10    nuv = normal_uncertain()
    11 
    12   where the required args are:
    13     descriptor    (str, description, '')
    14     mean          (float, mean, float('NaN'))
    15     stddev        (float, standard deviation, float('NaN'))
    16   and the optional args and defaults are:
    17     lower         (float, lower bound, -np.Inf)
    18     upper         (float, upper bound, np.Inf)
    19 
    20   note that zero arguments constructs a default instance, one
    21   argument of the class copies the instance, and three or more
    22   arguments constructs a new instance from the arguments.
    23 '''
    24     def __init__(self):
     12    """NORMAL_UNCERTAIN class definition
     13
     14    Usage:
     15        [nuv] = normal_uncertain(
     16            'descriptor', descriptor,
     17            'mean', mean,
     18            'stddev', stddev,
     19            'partition', partition
     20            )
     21
     22        where nuv is the normal_uncertain object returned by the constructor,
     23        mean and stddev are self explanatory, and partition is the partition
     24        vector for distributed variables. Can be a partition vector over
     25        elements or vertices.
     26
     27    Example:
     28        md.qmu.variables.rheology=normal_uncertain(
     29            'descriptor','RheologyBBar',
     30            'mean',1,
     31            'stddev',.05
     32            )
     33        md.qmu.variables.rheology=normal_uncertain(
     34            'descriptor','scaled_RheologyBBar',
     35            'mean',1,
     36            'stddev',.05,
     37            'partition',vpartition
     38            )
     39    """
     40
     41    def __init__(self): #{{{
    2542        self.descriptor = ''
    26         self.mean = float('NaN')
    27         self.stddev = float('NaN')
    28         self.lower = -np.Inf
    29         self.upper = np.Inf
    30 
    31     @staticmethod
    32     def normal_uncertain(*args):
     43        self.mean       = np.nan
     44        self.stddev     = np.nan
     45        self.partition  = []
     46        self.nsteps     = 0
     47    #}}}
     48
     49    @staticmethod
     50    def normal_uncertain(*args): #{{{
    3351        nargin = len(args)
    3452
    … …  
    4260                nuv = args[0]
    4361            else:
    44                 raise RuntimeError('Object ' + str(args[0]) + ' is a ' + str(type(args[0])) + ' class object, not "normal_uncertain".')
    45 
    46         # not enough arguments
    47         elif nargin == 2:
    48             raise RuntimeError('Construction of "normal_uncertain" class object requires at least 3 inputs.')
    49 
    50     # create the object from the input
    51         else:
    52             # lines differ here in other classes / tests; see asizec problem in notes
     62                raise Exception('Object ' + str(args[0]) + ' is a ' + str(type(args[0])) + ' class object, not "normal_uncertain".')
     63
     64        # create the object from the input
     65        else:
    5366            nuv = normal_uncertain()
    54             nuv.descriptor = str(args[0])
    55             nuv.mean = args[1]
    56             nuv.stddev = args[2]
    57             if nargin >= 4:
    58                 nuv.lower = args[3]
    59             if nargin >= 5:
    60                 nuv.upper = args[4]
    61             if nargin > 5:
    62                 print('WARNING: normal_uncertain:extra_arg: Extra arguments for object of class ' + str(type(nuv)) + '.')
    63 
    64         return [nuv]
    65 
    66     def __repr__(self):
    67         # display an individual object
    68         string = '\n'
    69         string += 'class "normal_uncertain" object = \n'
    70         string += '    descriptor: ' + str(self.descriptor) + '\n'
    71         string += '          mean: ' + str(self.mean) + '\n'
    72         string += '        stddev: ' + str(self.stddev) + '\n'
    73         string += '         lower: ' + str(self.lower) + '\n'
    74         string += '         upper: ' + str(self.upper) + '\n'
     67
     68            #recover options:
     69            options = pairoptions(*args)
     70
     71            #initialize fields:
     72            nuv.descriptor = options.getfieldvalue('descriptor')
     73            nuv.mean       = options.getfieldvalue('mean')
     74            nuv.stddev     = options.getfieldvalue('stddev')
     75
     76            #if the variable is scaled, a partition vector should have been
     77            #supplied, and that partition vector should have as many partitions
     78            #as the mean and stddev vectors:
     79            if nuv.isscaled():
     80                nuv.partition = options.getfieldvalue('partition')
     81                nuv.nsteps = options.getfieldvalue('nsteps', 1)
     82                npart = qmupart2npart(nuv.partition)
     83                if npart != nuv.mean.shape[0]:
     84                    raise Exception("normal_uncertain constructor: for the scaled variable %s the row size of the mean field should be identical to the number of partitions" % nuv.descriptor)
     85                if npart != nuv.stddev.shape[0]:
     86                    raise Exception("normal_uncertain constructor: for the scaled variable %s the row size of the stddev field should be identical to the number of partitions" % nuv.descriptor)
     87                if nuv.nsteps != nuv.mean.shape[1]:
     88                    raise Exception("normal_uncertain constructor: for the scaled variable %s the col size of the mean field should be identical to the number of time steps" % nuv.descriptor)
     89                if nuv.nsteps != nuv.stddev.shape[1]:
     90                    raise Exception("normal_uncertain constructor: for the scaled variable %s the col size of the stddev field should be identical to the number of time steps" % nuv.descriptor)
     91
     92        return [nuv] # Always return a list, so we have something akin to a MATLAB single row matrix
     93    #}}}
     94
     95    def __repr__(self): #{{{
     96        string = '   normal uncertain variable: '
     97        string = "%s\n%s" % (string, fielddisplay(self, 'descriptor', 'name tag'))
     98        string = "%s\n%s" % (string, fielddisplay(self, 'mean', 'pdf mean'))
     99        string = "%s\n%s" % (string, fielddisplay(self, 'stddev', 'pdf standard deviation'))
     100        if self.partition != []:
     101            string = "%s\n%s" % (string, fielddisplay(self, 'partition', 'partition vector defining where sampling will occur'))
     102        string = "%s\n%s" % (string, fielddisplay(self, 'nsteps', 'number of time steps'))
    75103
    76104        return string
    77 
    78     # from here on, nuv is either a single, or a 1d vector of, normal_uncertain
    79 
    80     @staticmethod
    81     def prop_desc(nuv, dstr):
    82         if type(nuv) not in [list, np.ndarray]:
    83             if nuv.descriptor != '' or type(nuv.descriptor) != str:
    84                 desc = str(nuv.descriptor)
    85             elif dstr != '':
    86                 desc = str(dstr)
    87             else:
    88                 desc = 'nuv'
    89             return desc
    90 
     105    #}}}
     106
     107    def __len__(self): #{{{
     108        if type(self.mean) in [list, np.ndarray]:
     109            return len(self.mean)
     110        else:
     111            return 1
     112    #}}}
     113
     114    def checkconsistency(self, md, solution, analyses): #{{{
     115        md = checkfield(md, 'field', self.mean, 'fieldname', 'normal_uncertain.mean', 'NaN', 1, 'Inf', 1, '>=', 0)
     116        md = checkfield(md, 'field', self.stddev, 'fieldname', 'normal_uncertain.stddev', 'NaN', 1, 'Inf', 1, '>=', 0)
     117        if self.isscaled():
     118            if self.partition == []:
     119                raise Exception("normal_uncertain is a scaled variable, but it's missing a partition vector")
     120            #better have a partition vector that has as many partitions as stddev's size:
     121            if self.stddev.shape[0] != partition_npart(self.partititon):
     122                raise Exception("normal_uncertain error message: row size of stddev and partition size should be identical")
     123            if self.mean.shape[0] != partition_npart(self.partition):
     124                raise Exception("normal_uncertain error message: row size of mean and partition size should be identical")
     125            #we need as many steps in stddev and mean as there are in time steps
     126            if self.stddev.shape[1] != self.nsteps:
     127                raise Exception("normal_uncertain error message: col size of stddev and partition size should be identical")
     128            if self.mean.shape[1] != self.nsteps:
     129                raise Exception("normal_uncertain error message: col size of mean and partition size should be identical")
     130            md = checkfield(md, 'field', self.partition, 'fieldname', 'normal_uncertain.partition', 'NaN', 1, 'Inf', 1, '>=', -1, 'numel', [md.mesh.numberofvertices, md.mesh.numberofvertices])
     131            if self.partition.shape[1] > 1:
     132                raise Exception("normal_uncertain error message: partition should be a column vector")
     133            partcheck = np.unique(self.partition)
     134            partmin = min(partcheck)
     135            partmax = max(partcheck)
     136            if partmax < -1:
     137                raise Exception("normal_uncertain error message: partition vector's min value should be -1 (for no partition), or start at 0")
     138            nmax = max(md.mesh.numberofelements, md.mesh.numberofvertices)
     139            if partmax > nmax:
     140                raise Exception("normal_uncertain error message: partition vector's values cannot go over the number of vertices or elements")
     141    #}}}
     142
     143    #virtual functions needed by qmu processing algorithms
     144    #implemented:
     145
     146    @staticmethod
     147    def prop_desc(nuv, dstr): #{{{
    91148        desc = ['' for i in range(np.size(nuv))]
    92149        for i in range(np.size(nuv)):
    … …  
    101158
    102159        return desc
    103 
    104     @staticmethod
    105     def prop_initpt(nuv):
    106         initpt = []
    107         return initpt
    108 
    109     @staticmethod
    110     def prop_lower(nuv):
    111         if type(nuv) not in [list, np.ndarray]:
    112             return nuv.lower
    113 
    114         lower = np.zeros(np.size(nuv))
    115         for i in range(np.size(nuv)):
    116             lower[i] = nuv[i].lower
    117 
    118         lower = allequal(lower, -np.inf)
    119 
    120         return lower
    121 
    122     @staticmethod
    123     def prop_upper(nuv):
    124         if type(nuv) not in [list, np.ndarray]:
    125             return nuv.upper
    126 
    127         upper = np.zeros(np.size(nuv))
    128         for i in range(np.size(nuv)):
    129             upper[i] = nuv[i].upper
    130 
    131         upper = allequal(upper, -np.inf)
    132         return upper
    133 
    134     @staticmethod
    135     def prop_mean(nuv):
    136         if type(nuv) not in [list, np.ndarray]:
    137             return nuv.mean
    138 
     160    #}}}
     161
     162    @staticmethod
     163    def prop_mean(nuv): #{{{
    139164        mean = np.zeros(np.size(nuv))
    140165        for i in range(np.size(nuv)):
    141166            mean[i] = nuv[i].mean
    142 
    143167        return mean
    144 
    145     @staticmethod
    146     def prop_stddev(nuv):
    147         if type(nuv) not in [list, np.ndarray]:
    148             return nuv.stddev
    149 
     168    #}}}
     169
     170    @staticmethod
     171    def prop_stddev(nuv): #{{{
    150172        stddev = np.zeros(np.size(nuv))
    151173        for i in range(np.size(nuv)):
    152174            stddev[i] = nuv[i].stddev
    153 
    154175        return stddev
    155 
    156     @staticmethod
    157     def prop_initst(nuv):
    158         initst = []
    159         return initst
    160 
    161     @staticmethod
    162     def prop_stype(nuv):
     176    #}}}
     177
     178    @staticmethod
     179    def prop_lower(nuv): #{{{
     180        lower = []
     181        return lower
     182    #}}}
     183
     184    @staticmethod
     185    def prop_upper(nuv): #{{{
     186        upper = []
     187        return upper
     188    #}}}
     189
     190    #default
     191    @staticmethod
     192    def prop_abscissas(hbu): #{{{
     193        abscissas = []
     194        return abscissas
     195    #}}}
     196
     197    @staticmethod
     198    def prop_pairs_per_variable(hbu): #{{{
     199        pairs_per_variable = []
     200        return pairs_per_variable
     201    #}}}
     202
     203    @staticmethod
     204    def prop_counts(hbu): #{{{
     205        counts = []
     206        return counts
     207    #}}}
     208    @staticmethod
     209    def prop_initpt(nuv): #{{{
     210        initpt = []
     211        return initpt
     212    #}}}
     213
     214    @staticmethod
     215    def prop_initst(nuv): #{{{
     216        inist = []
     217        return inist
     218    #}}}
     219
     220    @staticmethod
     221    def prop_stype(nuv): #{{{
    163222        stype = []
    164223        return stype
    165 
    166     @staticmethod
    167     def prop_scale(nuv):
     224    #}}}
     225
     226    @staticmethod
     227    def prop_scale(nuv): #{{{
    168228        scale = []
    169229        return scale
     230    #}}}
     231
     232    #new methods:
     233    def isdistributed(self): #{{{
     234        if strncmp(self.descriptor, 'distributed_', 12):
     235            return True
     236        else:
     237            return False
     238    #}}}
     239   
     240    def isscaled(self): #{{{
     241        if strncmp(self.descriptor, 'scaled_', 7):
     242            return True
     243        else:
     244            return False
     245    #}}}
    170246
    171247    @staticmethod
    172248    def dakota_write(fidi, dvar):
     249        # possible namespace pollution, the above import seems not to work
     250        from vlist_write import vlist_write
    173251        # collect only the variables of the appropriate class
    174         nuv = [struc_class(i, 'normal_uncertain', 'nuv') for i in dvar]
    175 
    176     # possible namespace pollution, the above import seems not to work
    177         from vlist_write import vlist_write
    178     # write variables
    179         vlist_write(fidi, 'normal_uncertain', 'nuv', nuv)
     252        nuv = deepcopy(dvar)
     253        fields = fieldnames(nuv)
     254        for field in fields:
     255            if getattr(nuv, field)[0].__class__.__name__ != 'normal_uncertain':
     256                delattr(nuv, field)
     257        if len(nuv) > 0:
     258            vlist_write(fidi, 'normal_uncertain', 'nuv', nuv)
     259    #}}}
     260

  • issm/trunk/src/m/classes/qmu/response_function.m

    r24686 r25836  
     1%RESPONSE_FUNCTION class definition
    12%
    2 %  definition for the response_function class.
    3 %
    4 %  [rf]=response_function(varargin)
    5 %
    6 %  where the required varargin are:
    7 %    descriptor    (char, description, '')
    8 %  and the optional varargin and defaults are:
    9 %    respl         (double vector, response levels, [])
    10 %    probl         (double vector, probability levels, [])
    11 %    rell          (double vector, reliability levels, [])
    12 %    grell         (double vector, gen. reliability levels, [])
    13 %
    14 %  note that zero arguments constructs a default instance; one
    15 %  argument of the class copies the instance; and one or more
    16 %  arguments constructs a new instance from the arguments.
    17 %
    18 %  "Copyright 2009, by the California Institute of Technology.
    19 %  ALL RIGHTS RESERVED. United States Government Sponsorship
    20 %  acknowledged. Any commercial use must be negotiated with
    21 %  the Office of Technology Transfer at the California Institute
    22 %  of Technology.  (J. Schiermeier, NTR 47078)
    23 %
    24 %  This software may be subject to U.S. export control laws.
    25 %  By accepting this  software, the user agrees to comply with
    26 %  all applicable U.S. export laws and regulations. User has the
    27 %  responsibility to obtain export licenses, or other export
    28 %  authority as may be required before exporting such information
    29 %  to foreign countries or providing access to foreign persons."
    30 %
     3%   Usage:
     4%      rf=response_function('descriptor',descriptor,'response_levels',respl,...
     5%           'probability_levels',probl,'reliability_levels',rell,general_reliability_levels',grell,...
     6%           'partition',partition);
     7%      where rf is the response function object returned by the constructor. All options except the
     8%      descriptor are optional. A partition can be provided for scaled variables.
     9%
     10%   Example:
     11%      md.qmu.responses.maxvel=response_function('descriptor','MaxVel','response_levels',[0],...
     12%                              'probl',[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]');
     13
    3114classdef response_function
    3215    properties
    … …  
    3619        rell      =[];
    3720        grell     =[];
     21        partition = [];
    3822    end
     23    methods
     24                function self=response_function(varargin) %constructor {{{
    3925
    40     methods
    41         function [rf]=response_function(varargin)
     26                        %recover options:
     27                        options = pairoptions(varargin{:});
    4228
    43             switch nargin
    44 
    45                                         %  create a default object
    46                 case 0
    47 
    48                                                  %  copy the object or create the object from the input
    49                 otherwise
    50                     if  (nargin == 1) && isa(varargin{1},'response_function')
    51                         rf=varargin{1};
    52                     else
    53                         asizec=num2cell(array_size(varargin{1:min(nargin,1)}));
    54                         rf(asizec{:})=response_function;
    55                         clear asizec
    56 
    57                         if ischar(varargin{1})
    58                             varargin{1}=cellstr(varargin{1});
    59                         end
    60                         for i=1:numel(rf)
    61                             if (numel(varargin{1}) > 1)
    62                                 rf(i).descriptor=varargin{1}{i};
    63                             else
    64                                 rf(i).descriptor=[char(varargin{1}) string_dim(rf,i,'vector')];
    65                             end
    66                         end
    67 
    68                         if (nargin >= 2)
    69                             for i=1:numel(rf)
    70                                 rf(i).respl     =varargin{2};
    71                             end
    72                             if (nargin >= 3)
    73                                 for i=1:numel(rf)
    74                                     rf(i).probl     =varargin{3};
    75                                 end
    76                                 if (nargin >= 4)
    77                                     for i=1:numel(rf)
    78                                         rf(i).rell      =varargin{4};
    79                                     end
    80                                     if (nargin >= 5)
    81                                         for i=1:numel(rf)
    82                                             rf(i).grell     =varargin{5};
    83                                         end
    84 
    85                                         if (nargin > 5)
    86                                             warning('response_function:extra_arg',...
    87                                                 'Extra arguments for object of class ''%s''.',...
    88                                                 class(rf));
    89                                         end
    90                                     end
    91                                 end
    92                             end
    93                         end
    94                     end
    95             end
    96 
    97         end
    98 
    99         function []=disp(rf)
    100 
    101         %  display the object
    102 
     29                        %initialize fields:
     30                        self.descriptor=getfieldvalue(options,'descriptor');
     31                        self.respl=getfieldvalue(options,'response_levels',[]);
     32                        self.probl=getfieldvalue(options,'probability_levels',...
     33                                [0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
     34                        self.rell=getfieldvalue(options,'reliability_levels',[]);
     35                        self.grell=getfieldvalue(options,'general_reliability_levels',[]);
     36                       
     37                        %if the response is scaled, a partition vector should have been supplied.
     38                        if self.isscaled(),
     39                                self.partition=getfieldvalue(options,'partition');
     40                                npart=partition_npart(self.partition);
     41                        end
     42                end %}}}
     43        function []=disp(rf)% {{{
     44            %TODO: Convert the following to fielddisplay
     45            % display the object
    10346            disp(sprintf('\n'));
    10447            for i=1:numel(rf)
    … …  
    10952                disp(sprintf('         probl: %s'      ,string_vec(rf(i).probl)));
    11053                disp(sprintf('          rell: %s'      ,string_vec(rf(i).rell)));
    111                 disp(sprintf('         grell: %s\n'    ,string_vec(rf(i).grell)));
     54                disp(sprintf('         grell: %s'    ,string_vec(rf(i).grell)));
     55                if ~isempty(rf.partition),
     56                                        fielddisplay(rf,'partition','partition vector defining where response will be computed');
     57                end
     58                disp(sprintf('\n'));
    11259            end
    11360
    114         end
     61        end % }}}
     62        function [desc]  =prop_desc(rf,dstr)% {{{
    11563
    116         function [desc]  =prop_desc(rf,dstr)
    11764            desc=cell(1,numel(rf));
    11865            for i=1:numel(rf)
    … …  
    12875            end
    12976            desc=allempty(desc);
    130         end
    131         function [stype] =prop_stype(rf)
     77        end % }}}
     78        function [stype] =prop_stype(rf)% {{{
     79
    13280            stype={};
    133         end
    134         function [scale] =prop_scale(rf)
     81        end% }}}
     82        function [scale] =prop_scale(rf)% {{{
     83
    13584            scale=[];
    136         end
    137         function [weight]=prop_weight(rf)
     85        end% }}}
     86        function [weight]=prop_weight(rf) % {{{
     87
    13888            weight=[];
    139         end
    140         function [lower] =prop_lower(rf)
     89        end% }}}
     90        function [lower] =prop_lower(rf)% {{{
     91
    14192            lower=[];
    142         end
    143         function [upper] =prop_upper(rf)
     93        end% }}}
     94        function [upper] =prop_upper(rf)% {{{
     95
    14496            upper=[];
    145         end
    146         function [target]=prop_target(rf)
     97        end % }}}
     98        function [target]=prop_target(rf) % {{{
    14799            target=[];
    148         end
    149         function [respl,probl,rell,grell]=prop_levels(rf)
     100        end % }}}
     101        function [respl,probl,rell,grell]=prop_levels(rf) % {{{
    150102            respl=cell(1,numel(rf));
    151103            probl=cell(1,numel(rf));
    … …  
    162114            rell =allempty(rell);
    163115            grell=allempty(grell);
    164         end
     116        end % }}}
     117                %new methods:
     118        function scaled =isscaled(self) % {{{
     119                if strncmp(self.descriptor,'scaled_',7),
     120                        scaled=1;
     121                else
     122                        scaled=0;
     123                end
     124        end % }}}
    165125    end
     126    methods (Static)
     127        function [rdesc]=dakota_write(fidi,dresp,rdesc) % {{{
    166128
    167     methods (Static)
    168         function [rdesc]=dakota_write(fidi,dresp,rdesc)
    169 
    170                           %  collect only the responses of the appropriate class
     129                        %collect only the responses of the appropriate class
    171130            rf=struc_class(dresp,'response_function');
    172 
    173                                 %  write responses
    174                                 [rdesc]=rlist_write(fidi,'response_functions','response_function',rf,rdesc);
    175                         end
    176 
    177                         function []=dakota_rlev_write(fidi,dresp,params)
     131                        %write responses
     132                        [rdesc]=rlist_write(fidi,'response_functions','response_function',rf,rdesc);
     133                end % }}}
     134                function []=dakota_rlev_write(fidi,dresp,params) % {{{
    178135
    179136                                %  collect only the responses of the appropriate class
    … …  
    182139                                %  write response levels
    183140                                rlev_write(fidi,rf,params);
    184         end
     141        end % }}}
    185142    end
    186143end

  • issm/trunk/src/m/classes/qmu/response_function.py

    r24313 r25836  
    11import numpy as np
     2
     3from fielddisplay import fielddisplay
     4from MatlabFuncs import *
     5from pairoptions import pairoptions
     6from partition_npart import *
     7from rlev_write import *
    28#from rlist_write import *
    3 from rlev_write import *
    4 from MatlabArray import *
    5 #move this later
    6 from helpers import *
    79
    810
    911class response_function(object):
    1012    '''
    11   definition for the response_function class.
    12 
    13   [rf] = response_function.response_function(args)
    14    rf = response_function()
    15 
    16   where the required args are:
    17     descriptor    (char, description, '')
    18   and the optional args and defaults are:
    19     respl         (double vector, response levels, [])
    20     probl         (double vector, probability levels, [])
    21     rell          (double vector, reliability levels, [])
    22     grell         (double vector, gen. reliability levels, [])
    23 
    24   note that zero arguments constructs a default instance, one
    25   argument of the class copies the instance, and one or more
    26   arguments constructs a new instance from the arguments.
    27 '''
    28 
     13    RESPONSE_FUNCTION class definition
     14
     15        Usage:
     16            rf = response_function(
     17                'descriptor', descriptor,
     18                'response_levels', respl,
     19                'probability_levels', probl,
     20                'reliability_levels',rell,
     21                'general_reliability_levels', grell,
     22                'partition', partition
     23            )
     24
     25            where rf is the response function object returned by the
     26            constructor. All options except the descriptor are optional. A partition can be provided for scaled variables.
     27
     28        Example:
     29            md.qmu.responses.maxvel = response_function(
     30                'descriptor', 'MaxVel',
     31                'response_levels',[0],
     32                'probl', [0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]
     33                )
     34    '''
    2935    def __init__(self):
    3036        self.descriptor = ''
    31         self.respl = []
    32         self.probl = []
    33         self.rell = []
    34         self.grell = []
     37        self.respl      = []
     38        self.probl      = []
     39        self.rell       = []
     40        self.grell      = []
     41        self.partition  = []
    3542
    3643    @staticmethod
    … …  
    4350
    4451        # copy the object or create the object from the input
    45         else:
    46             if nargin == 1 and isinstance(args[0], response_function):
     52        elif nargin == 1:
     53            if isinstance(args[0], response_function):
    4754                rf = args[0]
    4855            else:
    49                 asizec = array_size(*args[0:min(nargin, 1)])
    50                 rf = [response_function() for i in range(asizec[0]) for j in range(asizec[1])]
    51 
    52                 for i in range(np.size(rf)):
    53                     if (np.size(args[0]) > 1):
    54                         rf[i].descriptor = args[0][i]
    55                     else:
    56                         rf[i].descriptor = str(args[0]) + string_dim(rf, i, 'vector')
    57 
    58                 if nargin >= 2:
    59                     for i in range(np.size(rf)):
    60                         rf[i].respl = args[1]
    61 
    62                 if nargin >= 3:
    63                     for i in range(np.size(rf)):
    64                         rf[i].probl = args[2]
    65 
    66                 if nargin >= 4:
    67                     for i in range(np.size(rf)):
    68                         rf[i].rell = args[3]
    69 
    70                 if nargin >= 5:
    71                     for i in range(np.size(rf)):
    72                         rf[i].grell = args[4]
    73 
    74                 if nargin > 5:
    75                     print('WARNING: response_function:extra_arg: Extra arguments for object of class ' + str(type(rf)) + '.')
    76 
    77         return rf
    78 
    79     def __repr__(self):
    80         #  display the object
    81         string = '\n'
    82         string += 'class "response_function" object = \n'
    83         string += '    descriptor: ' + str(self.descriptor) + '\n'
    84         string += '         respl: ' + str(self.respl) + '\n'
    85         string += '         probl: ' + str(self.probl) + '\n'
    86         string += '          rell: ' + str(self.rell) + '\n'
    87         string += '         grell: ' + str(self.grell) + '\n'
     56                raise RuntimeError('Object ' + str(args[0]) + ' is a ' + str(type(args[0])) + ' class object, not "response_function".')
     57        else:
     58            # asizec = array_size(*args[0:min(nargin, 1)])
     59            # rf = [response_function() for i in range(asizec[0]) for j in range(asizec[1])]
     60
     61            # for i in range(np.size(rf)):
     62            #     if (np.size(args[0]) > 1):
     63            #         rf[i].descriptor = args[0][i]
     64            #     else:
     65            #         rf[i].descriptor = str(args[0]) + string_dim(rf, i, 'vector')
     66
     67            # if nargin >= 2:
     68            #     for i in range(np.size(rf)):
     69            #         rf[i].respl = args[1]
     70
     71            # if nargin >= 3:
     72            #     for i in range(np.size(rf)):
     73            #         rf[i].probl = args[2]
     74
     75            # if nargin >= 4:
     76            #     for i in range(np.size(rf)):
     77            #         rf[i].rell = args[3]
     78
     79            # if nargin >= 5:
     80            #     for i in range(np.size(rf)):
     81            #         rf[i].grell = args[4]
     82
     83            # if nargin > 5:
     84            #     print('WARNING: response_function:extra_arg: Extra arguments for object of class ' + str(type(rf)) + '.')
     85
     86            rf = response_function()
     87
     88            #recover options:
     89            options = pairoptions(*args)
     90
     91            #initialize fields:
     92            rf.descriptor = options.getfieldvalue('descriptor')
     93            rf.respl = options.getfieldvalue('response_levels', [])
     94            rf.probl = options.getfieldvalue('probability_levels', [0.0001, 0.001, 0.01, 0.25, 0.5, 0.75, 0.99, 0.999, 0.9999])
     95            rf.rell = options.getfieldvalue('reliability_levels', [])
     96            rf.grell = options.getfieldvalue('general_reliability_levels', [])
     97
     98            #if the response is scaled, a partition vector should have been supplied.
     99            if rf.isscaled():
     100                rf.partition = options.getfieldvalue('partition')
     101                npart = partition_npart(rf.partition)
     102
     103        return [rf] # Always return a list, so we have something akin to a MATLAB single row matrix
     104
     105    def __repr__(rf): #{{{
     106        # display the object
     107        string = 'class "response_function" object = \n'
     108        string = "%s\n%s" % (string, fielddisplay(rf, 'descriptor', 'name tag'))
     109        string = "%s\n%s" % (string, fielddisplay(rf, 'respl', 'response levels'))
     110        string = "%s\n%s" % (string, fielddisplay(rf, 'probl', 'probability levels'))
     111        string = "%s\n%s" % (string, fielddisplay(rf, 'rell', 'reliability levels'))
     112        string = "%s\n%s" % (string, fielddisplay(rf, 'grell', 'general reliability levels'))
     113
     114        if rf.partition != []:
     115            string = "%s\n%s" % (string, fielddisplay(rf, 'partition', 'partition vector defining where the response will be computed'))
    88116
    89117        return string
    90 
    91     def __len__(self):
     118    #}}}
     119
     120    def __len__(self): #{{{
    92121        return max(len(self.respl), len(self.probl), len(self.rell), len(self.grell))
    93 
    94     # from here on, rf is either a single, or a 1d vector of, response_function
    95 
    96     @staticmethod
    97     def prop_desc(rf, dstr):
    98         # response_function is always a vector, or should be, even with just 1
    99         if type(rf) not in [list, np.ndarray]:
    100             rf = [rf]
    101 
     122    #}}}
     123
     124    @staticmethod
     125    def prop_desc(rf, dstr): #{{{
    102126        desc = ['' for i in range(np.size(rf))]
    103127        for i in range(np.size(rf)):
    … …  
    111135        desc = allempty(desc)
    112136        return desc
    113 
    114     @staticmethod
    115     def prop_stype(rf):
     137    #}}}
     138
     139    @staticmethod
     140    def prop_stype(rf): #{{{
    116141        stype = []
    117142        return stype
    118 
    119     @staticmethod
    120     def prop_scale(rf):
     143    #}}}
     144
     145    @staticmethod
     146    def prop_scale(rf): #{{{
    121147        scale = []
    122148        return scale
    123 
    124     @staticmethod
    125     def prop_weight(rf):
     149    #}}}
     150
     151    @staticmethod
     152    def prop_weight(rf): #{{{
    126153        weight = []
    127154        return weight
    128 
    129     @staticmethod
    130     def prop_lower(rf):
     155    #}}}
     156
     157    @staticmethod
     158    def prop_lower(rf): #{{{
    131159        lower = []
    132160        return lower
    133 
    134     @staticmethod
    135     def prop_upper(rf):
     161    #}}}
     162
     163    @staticmethod
     164    def prop_upper(rf): #{{{
    136165        upper = []
    137166        return upper
    138 
    139     @staticmethod
    140     def prop_target(rf):
     167    #}}}
     168
     169    @staticmethod
     170    def prop_target(rf): #{{{
    141171        target = []
    142172        return target
     173    #}}}
    143174
    144175    @staticmethod
    145176    def prop_levels(rf):
    146         # response_function is always a vector, or should be, even with just 1
    147         if type(rf) not in [list, np.ndarray]:
    148             rf = [rf]
    149 
    150177        respl = empty_nd_list(np.size(rf))
    151178        probl = empty_nd_list(np.size(rf))
    152179        rell = empty_nd_list(np.size(rf))
    153180        grell = empty_nd_list(np.size(rf))
    154 
    155181        for i in range(np.size(rf)):
    156182            respl[i] = rf[i].respl
    … …  
    158184            rell[i] = rf[i].rell
    159185            grell[i] = rf[i].grell
    160 
    161186        respl = allempty(respl)
    162187        probl = allempty(probl)
    … …  
    165190        return [respl, probl, rell, grell]
    166191
     192    #new methods:
     193    def isscaled(self): #{{{
     194        if strncmpi(self.descriptor, 'scaled_', 7):
     195            return True
     196        else:
     197            return False
     198    #}}}
     199
    167200    @staticmethod
    168201    def dakota_write(fidi, dresp, rdesc):
    169         # collect only the responses of the appropriate class
    170         rf = [struc_class(vars(dresp)[i][j], 'response_function', 'rf') for i in fieldnames(dresp) for j in range(len(vars(dresp)[i]))]
    171 
    172202        #possible namespace pollution here
    173203        from rlist_write import rlist_write
    174         # write responses
    175         rdesc = rlist_write(fidi, 'response_function', 'rf', rf, rdesc)
    176 
     204        # collect only the responses of the appropriate class
     205        #rf = [struc_class(vars(dresp)[i][j], 'response_function', 'rf') for i in fieldnames(dresp) for j in range(len(vars(dresp)[i]))]
     206        resp = deepcopy(dresp)
     207        fields = fieldnames(resp)
     208        for field in fields:
     209            if getattr(resp, field)[0].__class__.__name__ != 'response_function':
     210                delattr(resp, field)
     211        if len(resp) > 0:
     212            rdesc = rlist_write(fidi, 'response_function', 'rf', resp, rdesc)
    177213        return rdesc
    178214

  • issm/trunk/src/m/classes/qmu/uniform_uncertain.m

    r24686 r25836  
     1%UNIFORM_UNCERTAIN Class definition
    12%
    2 %  definition for the uniform_uncertain class.
    3 %
    4 %  [uuv]=uniform_uncertain(varargin)
    5 %
    6 %  where the required varargin are:
    7 %    descriptor    (char, description, '')
    8 %    lower         (double, lower bound, -Inf)
    9 %    upper         (double, upper bound,  Inf)
    10 %
    11 %  note that zero arguments constructs a default instance; one
    12 %  argument of the class copies the instance; and three or more
    13 %  arguments constructs a new instance from the arguments.
    14 %
    15 %  "Copyright 2009, by the California Institute of Technology.
    16 %  ALL RIGHTS RESERVED. United States Government Sponsorship
    17 %  acknowledged. Any commercial use must be negotiated with
    18 %  the Office of Technology Transfer at the California Institute
    19 %  of Technology.  (J. Schiermeier, NTR 47078)
    20 %
    21 %  This software may be subject to U.S. export control laws.
    22 %  By accepting this  software, the user agrees to comply with
    23 %  all applicable U.S. export laws and regulations. User has the
    24 %  responsibility to obtain export licenses, or other export
    25 %  authority as may be required before exporting such information
    26 %  to foreign countries or providing access to foreign persons."
    27 %
     3%%   Usage:
     4%      nuv=uniform_uncertain('descriptor',descriptor,'lower',lower,'upper',upper,'partition',partition);
     5%      where nuv is the uniform_uncertain object returned by the constructor, lower and upper are the
     6%      pdf distribution bounds, and partition is the partition vector for distributed variables.
     7%      Can be a partition %      vector over elements or vertices.
     8%
     9%   Example:
     10%      md.qmu.variables.rheology=uniform_uncertain('descriptor','RheologyBBar','lower',1e8,'upper',1e9);
     11%      md.qmu.variables.rheology=uniform_uncertain('descriptor','RheologyBBar','lower',1e8,'upper',1e9,'partition',vpartition);
     12%
     13
    2814classdef uniform_uncertain
    29     properties
    30         descriptor='';
    31         lower     =-Inf;
    32         upper     = Inf;
    33     end
     15        properties
     16                descriptor      = '';
     17                lower           = -Inf;
     18                upper           = Inf;
     19                partition       = [];
     20                nsteps          = 0;
     21        end
     22        methods
     23                function self=uniform_uncertain(varargin) %constructor {{{
    3424
    35     methods
    36         function [uuv]=uniform_uncertain(varargin)
     25                        %recover options:
     26                        options=pairoptions(varargin{:});
    3727
    38             switch nargin
     28                        %initialize fields:
     29                        self.descriptor=getfieldvalue(options,'descriptor');
     30                        self.upper=getfieldvalue(options,'upper');
     31                        self.lower=getfieldvalue(options,'lower');
    3932
    40 %  create a default object
     33                        %if the variable is scaled, a partition vector should have been
     34                        %supplied, and that partition vector should have as many partitions
     35                        %as the lower and upper vectors:
     36                        if self.isscaled(),
     37                                self.partition=getfieldvalue(options,'partition');
     38                                self.nsteps=getfieldvalue(options,'nsteps',1);
     39                                npart=qmupart2npart(self.partition);
     40                                if npart~=size(self.upper,1),
     41                                        error(['uniform_uncertain constructor: for the scaled variable' self.descriptor ' the row size of the upper field should be identical to the number of partitions']);
     42                                end
     43                                if npart~=size(self.lower,1),
     44                                        error(['uniform_uncertain constructor: for the scaled variable' self.descriptor ' the row size of the lower field should be identical to the number of partitions']);
     45                                end
     46                                if self.nsteps~=size(self.upper,2),
     47                                        error(['uniform_uncertain constructor: for the scaled variable ' self.descriptor ' the col size of the upper field should be identical to the number of time steps']);
     48                                end
     49                                if self.nsteps~=size(self.lower,2),
     50                                        error(['uniform_uncertain constructor: for the scaled variable ' self.descriptor ' the col size of the lower field should be identical to the number of time steps']);
     51                                end
     52                        end
    4153
    42                 case 0
     54                end %}}}
     55                function disp(self) % {{{
    4356
    44 %  copy the object
     57                        disp(sprintf('   uniform uncertain variable: '));
     58                        fielddisplay(self,'descriptor','name tag');
     59                        fielddisplay(self,'lower','pdf lower bound');
     60                        fielddisplay(self,'upper','pdf upper bound');
     61                        if ~isempty(self.partition),
     62                                fielddisplay(self,'partition','partition vector defining where sampling will occur');
     63                        end
     64                        fielddisplay(self,'nsteps','number of time steps');
     65                end
     66                %}}}
     67                function md=checkconsistency(self,md,solution,analyses) % {{{
    4568
    46                 case 1
    47                     if isa(varargin{1},'uniform_uncertain')
    48                         uuv=varargin{1};
    49                     else
    50                         error('Object ''%s'' is a ''%s'' class object, not ''%s''.',...
    51                             inputname(1),class(varargin{1}),'uniform_uncertain');
    52                     end
     69                        md = checkfield(md,'field',self.upper,'fieldname','uniform_uncertain.upper','NaN',1,'Inf',1,'>',self.lower,'numel',length(self.lower));
     70                        md = checkfield(md,'field',self.lower,'fieldname','uniform_uncertain.lower','NaN',1,'Inf',1,'<',self.upper,'numel',length(self.upper));
     71                        if self.isscaled(),
     72                                if isempty(self.partition),
     73                                        error('uniform_uncertain is a scaled variable, but it''s missing a partition vector');
     74                                end
     75                                %better have a partition vector that has as many partitions as upper and lower's size:
     76                                if size(self.upper,1)~=partition_npart(self.partition),
     77                                        error('uniform_uncertain error message: row size of upper and partition size should be identical');
     78                                end
     79                                if size(self.lower,1)~=partition_npart(self.partition),
     80                                        error('uniform_uncertain error message: row size of lower and partition size should be identical');
     81                                end
     82                                %we need as steps in upper and lower as there are time steps:
     83                                if size(self.upper,2)~=self.nsteps,
     84                                        error('uniform_uncertain error message: col size of upper and number of time steps should be identical');
     85                                end
     86                                if size(self.lower,2)~=self.nsteps,
     87                                        error('uniform_uncertain error message: col size of lower and number of time steps should be identical');
     88                                end
    5389
    54 %  not enough arguments
    55 
    56                 case 2
    57                     error('Construction of ''%s'' class object requires at least %d inputs.',...
    58                         'uniform_uncertain',3)
    59 
    60 %  create the object from the input
    61 
    62                 otherwise
    63                     asizec=num2cell(array_size(varargin{1:min(nargin,3)}));
    64                     uuv(asizec{:})=uniform_uncertain;
    65                     clear asizec
    66 
    67                     if ischar(varargin{1})
    68                         varargin{1}=cellstr(varargin{1});
    69                     end
    70                     for i=1:numel(uuv)
    71                         if (numel(varargin{1}) > 1)
    72                             uuv(i).descriptor=varargin{1}{i};
    73                         else
    74                             uuv(i).descriptor=[char(varargin{1}) string_dim(uuv,i,'vector')];
    75                         end
    76                         if (numel(varargin{2}) > 1)
    77                             uuv(i).lower     =varargin{2}(i);
    78                         else
    79                             uuv(i).lower     =varargin{2};
    80                         end
    81                         if (numel(varargin{3}) > 1)
    82                             uuv(i).upper     =varargin{3}(i);
    83                         else
    84                             uuv(i).upper     =varargin{3};
    85                         end
    86                     end
    87             end
    88 
    89         end
    90 
    91         function []=disp(uuv)
    92 
    93 %  display the object
    94 
    95             disp(sprintf('\n'));
    96             for i=1:numel(uuv)
    97                 disp(sprintf('class ''%s'' object ''%s%s'' = \n',...
    98                     class(uuv),inputname(1),string_dim(uuv,i)));
    99                 disp(sprintf('    descriptor: ''%s'''  ,uuv(i).descriptor));
    100                 disp(sprintf('         lower: %g'      ,uuv(i).lower));
    101                 disp(sprintf('         upper: %g\n'    ,uuv(i).upper));
    102             end
    103 
    104         end
    105 
    106         function [desc]  =prop_desc(uuv,dstr)
    107             desc=cell(1,numel(uuv));
    108             for i=1:numel(uuv)
    109                 if ~isempty(uuv(i).descriptor)
    110                     desc(i)=cellstr(uuv(i).descriptor);
    111                 elseif ~isempty(inputname(1))
    112                     desc(i)=cellstr([inputname(1) string_dim(uuv,i,'vector')]);
    113                 elseif exist('dstr','var')
    114                     desc(i)=cellstr([dstr         string_dim(uuv,i,'vector')]);
    115                 else
    116                     desc(i)=cellstr(['uuv'        string_dim(uuv,i,'vector')]);
    117                 end
    118             end
    119             desc=allempty(desc);
    120         end
    121         function [initpt]=prop_initpt(uuv)
    122             initpt=[];
    123         end
    124         function [lower] =prop_lower(uuv)
    125             lower=zeros(1,numel(uuv));
    126             for i=1:numel(uuv)
    127                 lower(i)=uuv(i).lower;
    128             end
    129             lower=allequal(lower,-Inf);
    130         end
    131         function [upper] =prop_upper(uuv)
    132             upper=zeros(1,numel(uuv));
    133             for i=1:numel(uuv)
    134                 upper(i)=uuv(i).upper;
    135             end
    136             upper=allequal(upper, Inf);
    137         end
    138         function [mean]  =prop_mean(uuv)
    139             mean=[];
    140         end
    141         function [stddev]=prop_stddev(uuv)
    142             stddev=[];
    143         end
    144         function [initst]=prop_initst(uuv)
    145             initst=[];
    146         end
    147         function [stype] =prop_stype(uuv)
    148             stype={};
    149         end
    150         function [scale] =prop_scale(uuv)
    151             scale=[];
    152         end
    153                 function [abscissas] =prop_abscissas(hbu) % {{{
    154             abscissas=[];
    155         end % }}}
    156         function [counts] =prop_counts(hbu) % {{{
    157             counts=[];
    158         end % }}}
    159         function [pairs_per_variable] =prop_pairs_per_variable(hbu) % {{{
     90                                md = checkfield(md,'field',self.partition,'fieldname','uniform_uncertain.partition','NaN',1,'Inf',1,'>=',-1,'numel',[md.mesh.numberofvertices,md.mesh.numberofelements]);
     91                                if size(self.partition,2)>1,
     92                                        error('uniform_uncertain error message: partition should be a column vector');
     93                                end
     94                                partcheck=unique(self.partition);
     95                                partmin=min(partcheck);
     96                                partmax=max(partcheck);
     97                                if partmax<-1,
     98                                        error('uniform_uncertain error message: partition vector''s min value should be -1 (for no partition), or start at 0');
     99                                end
     100                                nmax=max(md.mesh.numberofelements,md.mesh.numberofvertices);
     101                                if partmax>nmax,
     102                                        error('uniform_uncertain error message: partition vector''s values cannot go over the number of vertices or elements');
     103                                end
     104                        end
     105                end % }}}
     106                %virtual functions needed by qmu processing algorithms:
     107                %implemented:
     108                function [desc]=prop_desc(uuv,dstr) % {{{
     109                        desc=cell(1,numel(uuv));
     110                        for i=1:numel(uuv)
     111                                if ~isempty(uuv(i).descriptor)
     112                                        desc(i)=cellstr(uuv(i).descriptor);
     113                                elseif ~isempty(inputname(1))
     114                                        desc(i)=cellstr([inputname(1) string_dim(uuv,i,'vector')]);
     115                                elseif exist('dstr','var')
     116                                        desc(i)=cellstr([dstr         string_dim(uuv,i,'vector')]);
     117                                else
     118                                        desc(i)=cellstr(['uuv'        string_dim(uuv,i,'vector')]);
     119                                end
     120                        end
     121                        desc=allempty(desc);
     122                end %}}}
     123                function [lower]=prop_lower(uuv) % {{{
     124                        lower=zeros(1,numel(uuv));
     125                        for i=1:numel(uuv)
     126                                lower(i)=uuv(i).lower;
     127                        end
     128                        lower=allequal(lower,-Inf);
     129                end %}}}
     130                function [upper]=prop_upper(uuv) % {{{
     131                        upper=zeros(1,numel(uuv));
     132                        for i=1:numel(uuv)
     133                                upper(i)=uuv(i).upper;
     134                        end
     135                        %upper=allequal(upper, Inf);
     136                end % }}}
     137                %default
     138                function [stddev]=prop_stddev(uuv)  %{{{
     139                        stddev=[];
     140                end % }}}
     141                function [mean]=prop_mean(nuv) % {{{
     142                        mean=[];
     143                end % }}}
     144                function [initpt]=prop_initpt(uuv) %{{{
     145                        initpt=[];
     146                end %}}}
     147                function [initst]=prop_initst(uuv) %{{{
     148                        initst=[];
     149                end %}}}
     150                function [stype]=prop_stype(uuv) %{{{
     151                        stype={};
     152                end %}}}
     153                function [scale]=prop_scale(uuv) %{{{
     154                        scale=[];
     155                end %}}}
     156                function [abscissas]=prop_abscissas(hbu) % {{{
     157                        abscissas=[];
     158                end % }}}
     159                function [counts]=prop_counts(hbu) % {{{
     160                        counts=[];
     161                end % }}}
     162                function [pairs_per_variable]=prop_pairs_per_variable(hbu) % {{{
    160163                        pairs_per_variable=[];
    161         end % }}}
    162 
    163     end
    164 
    165     methods (Static)
    166         function []=dakota_write(fidi,dvar)
    167 
    168 %  collect only the variables of the appropriate class
    169 
    170             uuv=struc_class(dvar,'uniform_uncertain');
    171 
    172 %  write variables
    173 
    174             vlist_write(fidi,'uniform_uncertain','uuv',uuv);
    175         end
    176     end
     164                end % }}}
     165                %new methods:
     166                        function distributed=isdistributed(self) % {{{
     167                        if strncmp(self.descriptor,'distributed_',12),
     168                                distributed=1;
     169                        else
     170                                distributed=0;
     171                        end
     172                end % }}}
     173        function scaled=isscaled(self) % {{{
     174                        if strncmp(self.descriptor,'scaled_',7),
     175                                scaled=1;
     176                        else
     177                                scaled=0;
     178                        end
     179                end % }}}
     180        end
     181        methods (Static)
     182                function []=dakota_write(fidi,dvar) % {{{
     183                        %  collect only the variables of the appropriate class
     184                        uuv=struc_class(dvar,'uniform_uncertain');
     185                        %  write variables
     186                        vlist_write(fidi,'uniform_uncertain','uuv',uuv);
     187                end %}}}
     188        end
    177189end

  • issm/trunk/src/m/classes/qmu/uniform_uncertain.py

    r24313 r25836  
    11import numpy as np
    2 #from vlist_write import *
     2
    33from MatlabArray import *
     4from MatlabFuncs import *
     5from fielddisplay import fielddisplay
     6from pairoptions import pairoptions
     7from partition_npart import *
     8from qmupart2npart import qmupart2npart
    49
    510
    611class uniform_uncertain(object):
    712    '''
    8   definition for the uniform_uncertain class.
    9 
    10   [uuv] = uniform_uncertain.uniform_uncertain(args)
    11    uuv = uniform_uncertain()
    12 
    13   where the required args are:
    14     descriptor    (str, description, '')
    15     lower         (float, lower bound, -np.Inf)
    16     upper         (float, upper bound, np.Inf)
    17 
    18   note that zero arguments constructs a default instance, one
    19   argument of the class copies the instance, and three or more
    20   arguments constructs a new instance from the arguments.
    21 '''
    22 
    23     def __init__(self):
     13    UNIFORM_UNCERTAIN class definition
     14
     15    Usage:
     16        [uuv] = uniform_uncertain(
     17            'descriptor', descriptor,
     18            'lower', lower,
     19            'upper', upper,
     20            'partition', partition
     21            )
     22
     23        where uuv is the uniform_uncertain object returned by the constructor,
     24        lower and upper are the pdf distribution bounds, and partition is the
     25        partition vector for distributed variables. Can be a partition vector
     26        over elements or vertices.
     27
     28    Example:
     29        md.qmu.variables.rheology = uniform_uncertain(
     30            'descriptor', 'RheologyBBar',
     31            'lower', 1e8,
     32            'upper', 1e9
     33            )
     34        md.qmu.variables.rheology = uniform_uncertain(
     35            'descriptor', 'RheologyBBar',
     36            'lower', 1e8,
     37            'upper', 1e9,
     38            'partition', vpartition
     39            )
     40    '''
     41    def __init__(self): #{{{
    2442        self.descriptor = ''
    25         self.lower = -np.Inf
    26         self.upper = np.Inf
    27 
    28     @staticmethod
    29     def uniform_uncertain(*args):
     43        self.lower      = -np.inf
     44        self.upper      = np.inf
     45        self.partition  = []
     46        self.nsteps     = 0
     47    #}}}
     48
     49    @staticmethod
     50    def uniform_uncertain(*args): #{{{
    3051        nargin = len(args)
    3152
    … …  
    3960                uuv = args[0]
    4061            else:
    41                 raise RuntimeError('Object ' + str(args[0]) + ' is a ' + str(type(args[0])) + ' class object, not "uniform_uncertain".')
    42 
    43         # not enough arguments
    44         elif nargin == 2:
    45             raise RuntimeError('Construction of "uniform_uncertain" class object requires at least 3 inputs.')
     62                raise Exception('Object ' + str(args[0]) + ' is a ' + str(type(args[0])) + ' class object, not "uniform_uncertain".')
    4663
    4764        # create the object from the input
    4865        else:
    49             # leaving this here in case it becomes important in the future
    50             #asizec = array_size(*args[0:min(nargin, 3)])
    51             #uuv = [uniform_uncertain() for i in range(asizec[0]) for j in range(asizec[1])]
    5266            uuv = uniform_uncertain()
    53             uuv.descriptor = str(args[0])
    54             uuv.lower = args[1]
    55             uuv.upper = args[2]
    56         if (nargin > 3):
    57             print('WARNING: uniform_uncertain:extra_arg: Extra arguments for object of class ' + type(uuv) + '.')
    58 
    59         return [uuv]
    60 
    61     def __repr__(self):
    62         # display an individual object
    63         string = '\n'
    64         string += 'class "uniform_uncertain" object = \n'
    65         string += '    descriptor: ' + str(self.descriptor) + '\n'
    66         string += '         lower: ' + str(self.lower) + '\n'
    67         string += '         upper: ' + str(self.upper) + '\n'
     67
     68            #recover options:
     69            options = pairoptions(*args)
     70
     71            #initialize fields:
     72            uuv.descriptor = options.getfieldvalue('descriptor')
     73            uuv.lower      = options.getfieldvalue('lower')
     74            uuv.upper      = options.getfieldvalue('upper')
     75
     76            #if the variable is scaled, a partition vector should have been
     77            #supplied, and  that partition vector should have as many
     78            #partitions as the lower and upper vectors:
     79            if uuv.isscaled():
     80                uuv.partition = options.getfieldvalue('partition')
     81                uuv.nsteps = options.getfieldvalue('nsteps', 1)
     82                npart = qmupart2npart(uuv.partition)
     83                if npart != uuv.upper.shape[0]:
     84                    raise Exception("uniform_uncertain constructor: for the scaled variable %s the upper field is not currently a vector of values for all the partitions described in the partition vector" % uuv.descriptor)
     85                if npart != uuv.lower.shape[0]:
     86                    raise Exception("uniform_uncertain constructor: for the scaled variable %s the lower field is not currently a vector of values for all the partitions described in the partition vector" % uuv.descriptor)
     87                if uuv.nsteps != uuv.upper.shape[1]:
     88                    raise Exception("uniform_uncertain constructor: for the scaled variable %s the col size of the upper field should be identical to the number of time steps" % uuv.descriptor)
     89                if uuv.nsteps != uuv.lower.shape[1]:
     90                    raise Exception("uniform_uncertain constructor: for the scaled variable %s the col size of the lower field should be identical to the number of time steps" % uuv.descriptor)
     91
     92        return [uuv] # Always return a list, so we have something akin to a MATLAB single row matrix
     93    #}}}
     94
     95    def __repr__(self): #{{{
     96        string = '   uniform uncertain variable: '
     97        string = "%s\n%s" % (string, fielddisplay(self, 'descriptor', 'name tag'))
     98        string = "%s\n%s" % (string, fielddisplay(self, 'lower', 'pdf lower bound'))
     99        string = "%s\n%s" % (string, fielddisplay(self, 'upper', 'pdf upper bound'))
     100        if self.partition != []:
     101            string = "%s\n%s" % (string, fielddisplay(self, 'partition', 'partition vector defining where sampling will occur'))
     102        string = "%s\n%s" % (string, fielddisplay(self, 'nsteps', 'number of time steps'))
    68103
    69104        return string
    70 
    71     # from here on, uuv is either a single, or a 1d vector of, uniform_uncertain
    72 
    73     @staticmethod
    74     def prop_desc(uuv, dstr):
    75         if type(uuv) not in [list, np.ndarray]:
    76             if uuv.descriptor != '' or type(uuv.descriptor) != str:
    77                 desc = str(uuv.descriptor)
    78             elif dstr != '':
    79                 desc = str(dstr)
    80             else:
    81                 desc = 'uuv'
    82             return desc
    83 
     105    #}}}
     106
     107    def __len__(self): #{{{
     108        if type(self.lower) in [list, np.ndarray]:
     109            return len(self.lower)
     110        else:
     111            return 1
     112    #}}}
     113   
     114    def checkconsistency(self, md, solution, analyses): #{{{
     115        md = checkfield(md, 'field', self.upper, 'fieldname', 'uniform_uncertain.upper', 'NaN', 1, 'Inf', 1, '>', self.lower, 'numel', len(self.lower))
     116        md = checkfield(md, 'field', self.lower, 'fieldname', 'uniform_uncertain.upper', 'NaN', 1, 'Inf', 1, '<', self.upper, 'numel', len(self.upper))
     117        if self.isscaled():
     118            if self.partition == []:
     119                raise Exception("uniform_uncertain is a scaled variable, but it's missing a partition vector")
     120            #better have a partition vector that has as many partitions as
     121            #upper and lower's size:
     122            if self.upper.shape[0] != partition_npart(self.partititon):
     123                raise Exception("uniform_uncertain error message: row size of upper and partition size should be identical")
     124            if self.lower.shape[0] != partition_npart(self.partition):
     125                raise Exception("uniform_uncertain error message: row size of lower and partition size should be identical")
     126            #we need as steps in upper and lower as there are time steps
     127            if self.stddev.shape[1] != self.nsteps:
     128                raise Exception("uniform_uncertain error message: col size of upper and partition size should be identical")
     129            if self.mean.shape[1] != self.nsteps:
     130                raise Exception("uniform_uncertain error message: col size of lower and partition size should be identical")
     131            md = checkfield(md, 'field', self.partition, 'fieldname', 'uniform_uncertain.partition', 'NaN', 1, 'Inf', 1, '>=', -1, 'numel', [md.mesh.numberofvertices, md.mesh.numberofvertices])
     132            if self.partition.shape[1] > 1:
     133                raise Exception("uniform_uncertain error message: partition should be a column vector")
     134            partcheck = np.unique(self.partition)
     135            partmin = min(partcheck)
     136            partmax = max(partcheck)
     137            if partmax < -1:
     138                raise Exception("uniform_uncertain error message: partition vector's min value should be -1 (for no partition), or start at 0")
     139            nmax = max(md.mesh.numberofelements, md.mesh.numberofvertices)
     140            if partmax > nmax:
     141                raise Exception("uniform_uncertain error message: partition vector's values cannot go over the number of vertices or elements")
     142    #}}}
     143
     144    #virtual functions needed by qmu processing algorithms:
     145    #implemented:
     146
     147    @staticmethod
     148    def prop_desc(uuv, dstr): #{{{
    84149        desc = ['' for i in range(np.size(uuv))]
    85150        for i in range(np.size(uuv)):
    … …  
    94159
    95160        return desc
    96 
    97     @staticmethod
    98     def prop_initpt(uuv):
    99         initpt = []
    100         return initpt
    101 
    102     @staticmethod
    103     def prop_lower(uuv):
    104         if type(uuv) not in [list, np.ndarray]:
    105             return uuv.lower
    106 
     161    #}}}
     162
     163    @staticmethod
     164    def prop_stddev(uuv): #{{{
     165        stddev = []
     166        return stddev
     167    #}}}
     168
     169    @staticmethod
     170    def prop_mean(uuv): #{{{
     171        mean = []
     172        return mean
     173    #}}}
     174
     175    @staticmethod
     176    def prop_lower(uuv): #{{{
    107177        lower = np.zeros(np.size(uuv))
    108178        for i in range(np.size(uuv)):
    109179            lower[i] = uuv[i].lower
    110180
    111         lower = allequal(lower, -np.Inf)
     181        lower = allequal(lower, -np.inf)
    112182
    113183        return lower
    114 
    115     @staticmethod
    116     def prop_upper(uuv):
    117         if type(uuv) not in [list, np.ndarray]:
    118             return uuv.upper
    119 
     184    #}}}
     185
     186    @staticmethod
     187    def prop_upper(uuv): #{{{
    120188        upper = np.zeros(np.size(uuv))
    121189        for i in range(np.size(uuv)):
    122190            upper[i] = uuv[i].upper
    123191
    124         upper = allequal(upper, np.Inf)
     192        #upper = allequal(upper, np.inf)
    125193
    126194        return upper
    127 
    128     @staticmethod
    129     def prop_mean(uuv):
    130         mean = []
    131         return mean
    132 
    133     @staticmethod
    134     def prop_stddev(uuv):
    135         stddev = []
    136         return stddev
    137 
    138     @staticmethod
    139     def prop_initst(uuv):
     195    #}}}
     196
     197    @staticmethod
     198    def prop_abscissas(hbu): #{{{
     199        abscissas = []
     200        return abscissas
     201    #}}}
     202
     203    @staticmethod
     204    def prop_pairs_per_variable(hbu): #{{{
     205        pairs_per_variable = []
     206        return pairs_per_variable
     207    #}}}
     208
     209    @staticmethod
     210    def prop_counts(hbu): #{{{
     211        counts = []
     212        return counts
     213    #}}}
     214
     215    @staticmethod
     216    def prop_initpt(uuv): #{{{
     217        initpt = []
     218        return initpt
     219    #}}}
     220
     221    @staticmethod
     222    def prop_initst(uuv): #{{{
    140223        initst = []
    141224        return initst
    142 
    143     @staticmethod
    144     def prop_stype(uuv):
     225    #}}}
     226
     227    @staticmethod
     228    def prop_stype(uuv): #{{{
    145229        stype = []
    146230        return stype
    147 
    148     @staticmethod
    149     def prop_scale(uuv):
     231    #}}}
     232
     233    @staticmethod
     234    def prop_scale(uuv): #{{{
    150235        scale = []
    151236        return scale
    152 
    153     @staticmethod
    154     def dakota_write(fidi, dvar):
    155         # collect only the variables of the appropriate class
    156         uuv = [struc_class(i, 'uniform_uncertain', 'uuv') for i in dvar]
     237    #}}}
     238
     239    #new methods:
     240    def isscaled(self): #{{{
     241        if strncmp(self.descriptor, 'scaled_', 7):
     242            return True
     243        else:
     244            return False
     245    #}}}
     246
     247    @staticmethod
     248    def dakota_write(fidi, dvar): #{{{
    157249        # possible namespace pollution, the above import seems not to work
    158250        from vlist_write import vlist_write
    159         # write variables
    160         vlist_write(fidi, 'uniform_uncertain', 'uuv', uuv)
     251        # collect only the variables of the appropriate class
     252        uuv = deepcopy(dvar)
     253        fields = fieldnames(uuv)
     254        for field in fields:
     255            if getattr(uuv, field)[0].__class__.__name__ != 'uniform_uncertain':
     256                delattr(uuv, field)
     257        if len(uuv) > 0:
     258            vlist_write(fidi, 'uniform_uncertain', 'uuv', uuv)
     259    #}}}
Note: See TracChangeset for help on using the changeset viewer.